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Abstract 


This report is a documentation of Version 1 of the Goddard Earth Observing System 
(GEOS) General Circulation Model (GCM) developed by NASA’s Data Assimilation 
Office and the Climate and Radiation Branch. The report is separated into three parts. 
Part I describes the overall features of the GCM, including numerical schemes for the 
hydrodynamics, sub-grid scale physical parameterizations, and boundary conditions. 
Part II provides a comprehensive description of the diagnostics available within the 
model, including the methodology for their selection and retrieval. Part III concludes the 
report with a User’s Guide, which includes the location of frozen versions of the GEOS 
system, instructions on how to create user-defined applications for running the GEOS 
system, and an overview of established production applications and output utilities 
available to the communitiy. 


ill 


PMMDtfM PAGE BLANK NOT FILMED 



Contents 


List of Figures viii 

I GEOS-1 GCM: A Descriptive Documention 1 

1 Introduction and Model Lineage 1 

1.1 Relevant Model Documentation 1 

2 Atmospheric Dynamics 2 

2.1 Horizontal and Vertical Discretization 4 

2.2 Time Integration Scheme 5 

2.3 Smoothing / Filling 7 

3 Atmospheric Physics 9 

3.1 Moist Convective Processes 9 

3.1.1 Sub-grid and Large-scale Convection 9 

3.1.2 Cloud Formation 11 

3.2 Radiation 12 

3.2.1 Shortwave Radiation . . . 13 

3.2.2 Longwave Radiation 13 

3.2.3 Cloud- Radiation Interaction 14 

3.3 Turbulence 15 

3.3.1 Planetary Boundary Layer . 16 

4 Boundary Conditions 16 


v 


J3LANK NOT FILMED 



4.1 Orography 16 

4.2 Land-Sea-Ice Surface Type 17 

4.3 Sea Surface Temperature 17 

4.4 Sea Surface Roughness 21 

4.5 Albedo 21 

4.6 Sea Ice 21 

4.7 Snow Cover 21 

4.8 Land Surface 21 

II GEOS-1 GCM: Diagnostics 23 

5 Introduction 23 

6 Diagnostic Overview 23 

7 GEOS GCM Diagnostic Menu 26 

8 Diagnostic Description 30 

III GEOS-1 GCM: User’s Guide 65 

9 Introduction 65 

10 GEOS Superstructure 65 

10.1 GEOS Applications 65 

10.2 Subroutines GCMDRV & RESTART 67 


vi 



11 GEOS Production 71 

11.1 Production Namelists 73 

11.2 Phoenix Format 77 

11.3 Geopotential Heights on Sigma Surfaces 81 

11.4 Interpolation from Sigma to Pressure Surfaces 83 

11.5 Sea-Level Pressure 88 

12 GEOS Unix Script 90 

12.1 Creating the GEOS GCM Namelist 90 

12.2 Creating the DIAGSIZE Subroutine 91 

12.3 Remote copy frozen GEOS Library 92 

12.4 Compile User Application and Output Routines 93 

12.5 Load and Run GEOS Simulation 94 

References 97 


vii 



List of Figures 


1 Sigma-level pressures used in the 20-level GEOS-1 GGM 6 

2 Shapiro filter response function used in the 2° x 2.5° GEOS-1 GCM 8 

3 Comparison between the Lanczos and mth - order Shapiro filter response func- 
tions for m — 2, 4, and 8 

4 GEOS-1 GCM Land/Water Mask and Topography used at 4° x 5° resolution. 19 

5 GEOS-1 GCM Land/Water Mask and Topography used at 2° x 2.5° resolution. 20 

6 Vertical placement and index notation for sigma levels in the GEOS GCM . 83 

7 Index notation for vertical interpolation of Geopotential Heights 84 



Part I 


GEOS-l GCM: A Descriptive 
Documention 

1 Introduction and Model Lineage 


The Goddard Earth Observing System-1 (GEOS-l) General Circulation Model (GCM) 
was developed by the Data Assimilation Office (DAO) at the Goddard Laboratory for 
Atmospheres (GLA), in collaboration with the Climate and Radiation Branch, for use in 
the system being developed to analyze EOS data. The GEOS Data Assimilation System 
(DAS), currently using the GEOS-l GCM in conjuction with an Optimal Interpolation 
(01) analysis scheme, has been used to produce a multi-year global atmospheric data set 
for climate research (Schubert et al., 1993). The GEOS-l GCM has also been used to 
produce multiple 10-year climate simulations as part of the DAO’s participation in the 
Atmospheric Model Intercomparison Project (AMIP) sponsored by the Program for Climate 
Model Diagnostics and Intercomparison (PCMDI) (see Gates, 1992). The GEOS-l DAS has 
been used operationally to provide scientific flight guidance during NASA’s participation in 
the Measurements for Assessing the Effects of Stratospheric Aircraft (MAESA) experiment. 

The earliest predecessor of the GE10S-1 GCM was developed in 1989 based on the “plug- 
compatible” concepts outlined in Kalnay et al. (1989), and subsequently improved in 1991 
(Fox- Rabinovitz, et al., 1991, Helfand et al., 1991). The plug-compatibility of the physi- 
cal parameterizations together with the plug-compatible “Dynamical Core” introduced by 
Suarez and Takacs (1994) facilitate the development and testing of new algorithms. To- 
gether the DAO and the Climate and Radiation Branch at GLA have produced a library 
of physical parameterizations and dynamical algorithms which can be utilized for various 
GCM applications. 


1.1 Relevant Model Documentation 

A description of the overall features of the GEOS Data Assimilation System (DAS) used by 
the DAO may be found in Schubert et al. (1993). A comprehensive documentation of the 
Aries/ GEOS Dynamical Core incorporating the horizontal and vertical discretization and 
finite-difference schemes used is given in Suarez and Takacs (1994). The Relaxed Arakawa- 
Schubert cumulus convective parameterization and the re-evaporation of falling rain are 
based upon the works of Moorthi and Suarez (1992) and Sud and Molod (1988), while the 
radiative processes are described by Harshvardhan et al. (1987). The planetary boundary 
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layer (PBL) and the upper level turbulence parameterization are based on the level 2.5 clo- 
sure model of Helfand and Labraga (1988) and Helfand et al. (1991). Additional postscript 
copies of this document may be obtained through anonymous ftp from hera.gsfc.nasa.gov. 
The path and filename of the document is: 

pub / gem /geos 1 .0_gcm.doc.ps 


2 Atmospheric Dynamics 


The momentum equations used in the GEOS-1 GCM are written in the “vector invariant” 
form, as in Sadourney (1975) and Arakawa and Lamb (1981), to facilitate the derivation 
of the energy and potential enstrophy conserving differencing scheme. The thermodynamic 
(potential temperature) and moisture (specific humidity) equations are written in flux form 
to facilitate potential temperature and moisture conservation. A complete description of 
the finite-difference scheme used can be found in Suarez and Takacs (1994). 


The GEOS GCM uses a o coordinate defined by 


a — 


p-pT 

y 

7 r 


(i) 


where p is the pressure, ir = p s — pr , p s is the surface pressure, and pj is a constant 
prescribed pressure at the top of the model atmosphere. With px — 0 this coordinate 
reduces to the conventional a coordinate proposed by Phillips (1957). 


With this vertical coordinate, the continuity equation becomes 
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dt 


-Vo- • (ttv) - 


d(Tcr) 
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(2) 


where v is the horizontal velocity vector. Integrating (2) and assuming a - 0 at p = p s and 
p = pj>, we obtain the forms used in the model: 


dir 

dt 
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and 
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The equation of state for an ideal gas is a — RT / p , where a is the specific density, T is the 
temperature, and R is the gas constant. The following alternative forms will be used below 
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where 0 = T/P is the potential temperature, P = (p/po) K , k = R/c p , c p is the specific heat 
at constant pressure, and po is a reference pressure. In obtaining the forms in (5) we have 
used Or- — k— and the relation 

op p 

For the time being virtual effects are neglected. 


The hydrostatic equation is 


d<& 

dp 


= —ck, 


where $ is the geopotential. Using (5) and (6), this can be written: 
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From (7) we obtain 


d$ 


— Cp0, 


dP ~ p 

which, following Arakawa and Suarez (1983), is the form used in the model. 
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( 8 ) 


The thermodynamic equation is written in flux form to facilitate the derivation of a B- 
conserving differencing scheme: 


d(ird) 

dt 
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d(irdB) ttQ 


( 9 ) 


where Q is the diabatic heating per unit mass. 


In addition to the equations of motion, the Aries/GEOS Dynamical Core computes tenden- 
cies for an arbitrary number of atmospheric contituents, such as water vapor and ozone. 
These are also written in flux form: 


M = _v„.(W*>)-^3^ + -r5W, (10) 

where q ^ is the specific mass of the kth. constituent, and is its source per unit mass 
of air. 


The momentum equation is written in “vector-invariant” form, as in Sadourney (1975) 
and Arakawa and Lamb (1981), to facilitate the derivation of an energy- and enstrophy- 
conserving differencing scheme: 
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where 


„ = w±« 

7 r 

is an “external” potential vorticity, / is the Coriolis parameter, k is the unit vector in the 
vertical, 

( = V ff xv 

is the vertical component of the vorticity along a surfaces, 

K = |(v-v) 

is the kinetic energy per unit mass, g is the acceleration of gravity, and T is the horizontal 
frictional stress. 


2.1 Horizontal and Vertical Discretization 

The GEOS-1 GCM is constructed in the horizontal using the staggered Arakawa C-grid 
and employs Version 1 of the Aries/GEOS Dynamical Core for the finite-differencing al- 
gorithm. The Aries/GEOS Dynamical Core is a plug- compatible dynamics module which 
is used at both the DAO in the GEOS GCM as well as at the Climate and Radiation 
Branch in the Aries GCM. The Aries GCM has been extensively used in many climate and 
coupled ocean/ atmosphere simulations (eg. Schubert et al. 1993, Higgins and Schubert, 
1993). Version 1 of the Aries/GEOS Dynamical Core is the second-order energy and poten- 
tial enstrophy conserving scheme of Sadourny described by Burridge and Haseler (1977), 
while Version 2, to be used in subsequent versions of the GEOS GCM, is the fourth-order 
formulation described in Suarez and Takacs, 1994. This scheme conserves total energy 
and potential enstrophy for the non-divergent component of the flow in the shallow water 
equations. 

The Aries/GEOS Dynamical Core uses a Lorenz or unstaggered vertical grid in generalized 
sigma coordinates. The vertical differencing scheme is that of Arakawa and Suarez (1983) 
which ensures that: 


• The pressure gradient force generates no circulation of vertically in- 
tegrated momentum along a contour of surface topography 

• The finite- difference analogues of the energy- conversion term have 
the same form in the kinetic energy and thermodynamic equations 

• The global mass integral of the potential temperature is conserved 
under adiabatic processes 

• The hydrostatic equation for the lowest thickness has a local form 

• The hydrostatic equation is exact for vertically isentropic atmo- 
spheres 
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• The pressure-gradient force is exact for three-dimensionally isen- 
tropic atmospheres 


The GEOS-1 AMIP simulation was run using a global 4° x 5° latitude/longitude grid in 
the horizontal, together with 20 sigma levels in the vertical. The GEOS-1 DAS 5-year 
re-analysis was run using 2° x 2.5° horizontal resolution. The vertical distribution of the 
sigma levels is such as to provide enhanced resolution in the planetary boundary layer and 
at the jet level. The model top pressure was prescribed to be 10 mb with a rigid lid. For 
a surface pressure of 1000 mb, the lowest sigma level is at 994 mb while the highest sigma 
level is at 19 mb. There are 5 sigma levels below 800 mb, and 7 sigma levels above 200 mb. 
The sigma-level pressure used in the 20-level GEOS-1 GCM are shown in Figure 1. 


2.2 Time Integration Scheme 


The GEOS GCM has the ability to use the Matsuno time integration scheme or the Leapfrog 
time integration scheme together with an Asselin (1972) time filter. The GEOS GCM multi- 
year simulations have all been run using the Leapfrog time scheme together with the Asselin 
averaging parameter equal to 0.05. The GEOS-DAS 5-year Re-analysis used the Matsuno 
time scheme. The GEOS GCM employs a somewhat unique method for incorporating 
adjustments due to diabatic processes (ie, moist convection, radiation and turbulence) and 
the analysis increments during an assimilation. At every model time step, all prognostic 
fields are updated due to both dynamical and sub-grid scale diabatic processes, shown 
schematically using the Leapfrog time scheme for an arbitrary prognostic field q as: 


where 
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(14) 


with S 0 defined as the solar contant, R a as the earth-sun distance in Astronomical Units, 
and cos(j) z as the cosine of the zenith angle. The Dynamics and Shapiro Filter time ten- 
dencies are updated every model time step. The diabatic time tendencies are computed 
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igure 1: Sigma- level pressures used m the zu-levej 






(updated) at a time step appropriate to the physical parameterizations, and are held con- 
stant between Physics calls. The time tendency for Moist Convection is updated every 10 
minutes, and for Turbulence every 30 minutes. The time tendency for Longwave Radiation 
is updated every 3 hours. Shortwave Radiation is updated once every three hours assuming 
a normalized incident solar radiation, and adjusted at every model time step by the true 
incident radiation. During GEOS-1 DAS assimilations, the Analysis Increment is updated 
every synoptic time period, or 6 hours. By gradually incorporating the diabatic adjust- 
ments during the model integration, shocks and dynamical imbalances are greatly reduced 
(cf. Bloom, et al. 1991). 

2.3 Smoothing / Filling 

An eighth- (sixteenth-) order Shapiro (1970) filter adjustment is computed for the winds, 
potential temperature, and specific humidity in conjunction with the 2° x 2.5° (4° x 5° ) 
GEOS GCM in order to globally damp small-scale dispersive waves. In order to reduce 
the dynamical imbalances generated by intermittent use of the full filter, a Shapiro filter 
tendency for any quantity q is defined as: 


(|) , . (15) 

V Ol / Shapiro Filter T 

where q F is the quantity after application of the full Shapiro filter, q is the unfiltered 
quantity, and r is an adjustable time scale. Thus, a fraction ~ of the full Shapiro filter is 
incorporated at each model timestep. In the GEOS-1 GCM, the time scale r is equal to 1.5 
hours, and has been chosen so as to remove the two-grid interval wave in about six hours. 
The Shapiro filter response function used in the 2° x 2.5° GEOS-1 GCM is shown in Figure 
2. The GEOS-1 GCM also uses a high-latitude Fourier filter to avoid linear instability due 
to violation of the CFL condition for the Lamb wave and internal gravity waves near the 
poles. This polar filter is only applied to the time-tendencies of the prognostic fields, ie. 
winds, potential temperature, specific humidity, and surface pressure. 

For the moisture equation, negative values computationally generated through advection 
are filled by “borrowing” from below while conserving the vertically integrated moisture, 
ie., 


(f^pqdz) =(f le ~ l pqdz) 
\Jle+l ) final \Jle + 1 ) initial 


Using the hydrostatic relation pSz 


} f inal 

^ — —~6a, we may write 
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8th— Order Shapiro Filter Response Function 



Figure 2: Shapiro filter response function used in the 2° x 2.5° GEOS-1 GCM. 
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Approximating equaton (17) by 
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(18) 


an expression for the updated moisture at level / is given by 


*«/<»»! = { Tqi + Tqi ~ l ~K^r 


- I Kqi -\— - — 
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(19) 


Assuming that irqi-i is initially negative, we require that its final value is set to zero. Thus 

KQl-l final 
final 

This process is repeated until the lowest level is reached. If the resulting 7 rq n i ay is negative, 
the mass-weighted specific humidity in the lowest model level is simply set to zero. 


= o, 


= vqi + nqi-i 


Act/-] 
A a i 


initial 
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3 Atmospheric Physics 

3.1 Moist Convective Processes 

3.1.1 Sub-grid and Large-scale Convection 

Sub-grid scale penetrative and shallow cumulus convection are parameterized using the Re- 
laxed Arakawa Schubert (RAS) scheme of Moorthi and Suarez (1992), which is a simplified 
Arakawa Schubert type scheme. The RAS scheme predicts mass fluxes from cloud types 
which have different entrainment rates and levels of neutral bouyancy, depending on the 
properties of the large scale environment. 

The thermodynamic variables that are used in RAS to describe the grid scale vertical profile 
are the dry static energy, s — c p T + gz , the moist static energy, h = c p T -f gz + Lq , and the 
saturation moist static energy, h* = c p T + gz + Lq*, where the superscript * refers to the 
saturation value. 
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The conceptual model behind RAS depicts the cloud as a rising plume, entraining mass 
from the environment during ascent, and detraining mass at the level of neutral buoyancy. 
RAS assumes that the cloud mass flux, ?/, normalized by the cloud base mass flux, is a 
linear function of height, expressed as: 


drf(z) dr](P K ) 

dz dP K 



where we have used the hydrostatic equation written in the form: 


_ _ c JLa 

dP K q 


The entrainment parameter, A, characterizes a particular cloud type based on its detrain- 
ment level, and is obtained by assuming that the level of detrainment is the level of neutral 
buoyancy, ie., the level at which the moist static energy of the cloud, h c , is equal to the sat- 
uration moist static energy of the environment, h * . Following Moorthi and Suarez (1992), 
A may be written as 

A = Hb ~ h * D 

~ ffpSWi-w' 

where the subscript B refers to cloud base, and the subscript D refers to the detrainment 
level. 


The convective instability is measured in terms of the cloud work function A, defined as the 
rate of change of cumulus kinetic energy. The cloud work function is related to the buoyancy, 
or the difference between the moist static energy in the cloud and in the environment: 


A 




Pb 


p d 1 + 7 L P K 


hr-fl* 


dP« 


where 7 is obtained from the Claussius Clapeyron equation, and the subscript c refers 

to the value inside the cloud. 


To determine the cloud base mass flux, the rate of change of A in time due to dissipation by 
the clouds is assumed to approximately balance the rate of change of A due to the generation 
by the large scale. This is the quasi-equilibrium assumption, and results in an expression 
for ms- 


m B 


_ dA 


dt 

Is 


where K is the cloud kernel, defined as the rate of change of the cloud work function 
per unit cloud base mass flux, and is obtained by differentiating the expression for A in 
time. The rate of change of A due to the generation by the large scale is written as the 
difference between the current A(t T A t) and its value after the previous convective time 
step A(t), divided by the time step. Since the convective parameterization is designed to 
nearly neutralize the instability, A(t) is approximated as some (near zero) critical A cr i t . 
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The predicted convective mass fluxes are used to solve grid-scale temperature and mois- 
ture budget equations to determine the impact of convection on the large scale fields of 
temperature (through latent heating and compensating subsidence) and moisture (through 
precipitation and detrainment): 

30 i m B ds 

tt7 = 

dt \ c c p P K dp 


and 


dq 

dt 


m B ( dh ds 
= aR -r 9n{ dt ~ dv ] 


where 9 — 


fz, and P 


(. P/Po )• 


As an approximation to a full interaction between the different allowable cloud types, many 
clouds are simulated frequently, each modifying the large scale environment some fraction a 
of the total adjustment. The parameterization thereby “relaxes” the large scale environment 
towards neutrality. 


In addition to the RAS cumulus convection scheme, the GEOS-1 GCM employs a Kessler- 
type scheme for the re-evaporation of falling rain (Sud and Molod, 1988), which corre- 
spondingly adjusts the impact on the large scale environment by the factor R. The scheme 
accounts for the rainfall intensity, the drop size distribution, and the temperature, pressure 
and relative humidity of the surrounding air. The moisture deficit in any model layer into 
which the rain may re-evaporate is a free parameter. 

Due to the increased vertical resolution in the Planetary Boundary Layer (PBL), the lowest 
two model layers are averaged to provide the sub-cloud layer for RAS (nominally 50 mb 
thick). Each time RAS is invoked (every ten simulated minutes), the possiblity for shallow 
convection is checked for the two layers just above cloud base. RAS also randomly chooses 
10 other cloud-top levels for the possibility of convection, from just above cloud base to the 
model top layer. 


Supersaturation or large-scale convection is defined in the GEOS GCM whenever the specific 
humidity in any grid-box exceeds its supersaturation value. The large-scale precipitation 
scheme rains at supersaturation, and re-evaporates during descent to partially saturate 
lower layers in a process that accounts for some simple microphysics. 


3.1.2 Cloud Formation 

Convective and large-scale cloudiness which is used for cloud- radiative interactions are de- 
termined diagnostically as part of the cumulus and large-scale parameterizations. The 
convective and large-scale cloud fractions are combined into two separate arrays for use in 
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the shortwave and longwave radiation packages, one array for random overlap (CLRO) and 
one array for maximum overlap (CLMO) cloudiness. 

If convection occurs with a cloud-top pressure less than 400 mb, a cloud fraction equal to 
1.0 is assigned into CLMO at the detrainment level. These tall convective towers are said 
to be maximumly overlapped, ie. they are radiatively correlated in the vertical. If shallow 
convection occurs with a detrainment level within 1 or 2 levels above the cloud base, a cloud 
fraction equal to 0.5 is assigned into CLRO, again at the detrainment level. These shallow 
clouds are not correlated in the vertical and are said to be randomly overlapped. No cloud 
fractions are prescribed to convection whose cloud-top is between 400 mb and 2 levels above 
cloud-base. 

Supersaturation or large scale cloudiness is defined whenever the large scale precipitation 
scheme determines that the grid box at any level becomes supersaturated. In order to ensure 
that at any instant the total cloud fraction be less than or equal to 1.0, supersaturation 
clouds are only prescribed when there are no deep convective clouds. Under such conditions, 
the grid box is assumed to be instantaneously fully covered with randomly overlapped 
clouds, with a cloud fraction of 1.0 being assigned into CLRO. 


3.2 Radiation 

The parameterization of radiation includes both shortwave radiation and longwave radiation 
in the GEOS-1 GCM. A single model “sponge” layer is inserted into the GEOS-1 GCM above 
Pjop with a temperature equal to the GCM’s first level temperature. This is done in order 
to avoid exposing the first GCM level to a zero downward thermal flux condition at Pt op • 
Radiative fluxes are calculated at each model edge-level in both up and down directions. 
The heating rates/cooling rates are then obtained from the vertical divergence of the net 
radiative fluxes. 

The net flux is 

p — p\ _ pi 

where F is the net flux, F 1 is the upward flux and F 1 is the downward flux. 

The heating rate due to the divergence of the radiative flux is given by 

dpc p T _ dF_ 
dt dz 

or 

dT _ g OF 
dt c p 7T d<J 

where g is the accelation due to gravity and c p is the heat capacity of air at constant 
pressure. 
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The infrared and solar radiation parameterizations employed in the GEOS-1 GCM follow 
closely those described by Harshvardhan et al. (1987). The time tendency for Longwave 
Radiation is updated every 3 hours. The time tendency for Shortwave Radiation is updated 
once every three hours assuming a normalized incident solar radiation, and subsequently 
modified at every model time step by the true incident radiation. For the AMIP simulations, 
the GEOS-1 GCM used the AMIP-prescribed solar constant value of 1365 W/m 2 and a C0 2 
mixing ratio of 345 ppm. For the 5-year GEOS-DAS re-analysis, a solar constant value of 
1380 W/m 2 was used together with a C0 2 mixing ratio value of 330 ppm. For the ozone 
mixing ratio, monthly mean zonally averaged climatological values specified as a function of 
latitude and height (Rosenfield, et al., 1987) are linearly interpolated to the current time. 


3.2.1 Shortwave Radiation 


The parameterization of heating due to shortwave radiation is an extension of that developed 
by Lads and Hansen (1974). The documentation of the computational aspects of the scheme 
are given in Davis (1982). Harshvardhan improved the treatment of the cloud scattering 
by allowing it to depend on the solar zenith angle. Two major absorbers, H 2 0 and 0 3 , 
are considered in the scheme. Radiative heating by water vapor absorption is modeled in 
the 0.7-0. 4 fim region. The total absoption of the ozone includes the absoption of visible 
radiation by the Chappuis bands (450-750 nm) and the absorption of ultraviolet radiation 
by the Hartly (240-280 nm) and Huggins bands (280-360 nm). For a cloudy atmosphere, the 
Meadow- Weaver (1980) modified Eddington approximation is adopted to solve the radiative 
transfer equation in a scattering medium. The multiple-scattering computation with the 
k-distribution functions is used to calculate the cloudy-sky water absorption. 


3.2.2 Longwave Radiation 


The parameterization of longwave radiation employs a wide band model. Four broadband 
transmissions have been used in the model. The parameterization of the H 2 0 , containing 
the two water bands at 15 (im and 9.6 /im, is calculated using the method of Chou (1984) 
based on the water vapor line and water vapor e-type absorption approximations. The C0 2 
band employs the scheme of Chou and Peng (1983) which seperates the band wing and 
band center scaled paths. The O 3 parameterization approach applies the formulation of 
Rosenfield et al. (1987), which modifies the line width used by Rodgers (1968) to include 
the affects of Doppler broadening. 
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3.2.3 Cloud-Radiation Interaction 


The cloud fraction values produced by the Moist Convective Processes are used for cloud- 
radiation interactions as a means for evaluating clear line-of-site probabilities and effective 
optical thicknesses. 

If we define the time-averaged random and maximum overlapped cloudiness as CLRO and 
CLMO respectively, then the probability of clear sky associated with random overlapped 
clouds at any level is (1-CLRO) while the probability of clear sky associated with maximum 
overlapped clouds at any level is (1-CLMO). The total clear sky probability is given by 
(1-CLR0)*(1-CLM0), thus the total cloud fraction at each level may be obtained by 1(1- 
CLR0)*(1-CLM0). 

At any given level, we may define the clear line-of-site probability by appropriately account- 
ing for the maximum and random overlap cloudiness. The clear line-of-site probability is 
defined to be equal to the product of the clear line-of-site probabilities associated with 
random and maximum overlap cloudiness. The clear line-of-site probability C(p,p') associ- 
ated with maximum overlap clouds, from the current pressure p to the model top pressure, 
p ' = ptop, or the model surface pressure, p 1 — p SU rfi is simply 1.0 minus the largest maximum 
overlap cloud value along the line-of-site, ie. 


1 - MAX X ( CLMOp ) 

Thus, even in the time-averaged sense it is assumed that the maximum overlap clouds are 
correlated in the vertical. The clear line-of-site probability associated with random overlap 
clouds is defined to be the product of the clear sky probabilities at each level along the 
line-of-site, ie. 


p' 

H(1-CLRO p ) 

p 

The total cloud fraction at a given level associated with a line- of-site calculation is given 

by 


p 

1 - (l - MAX$' [C LM Op] J J] (1 - CLROp) 

P 

The GEOS-1 GCM cloud optical thicknesses are specified based on cloud type and tem- 
perature. The “maximum overlap” clouds are assigned an optical thickness of 0.16 km~ x 
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for the shortwave radiative feedback. The “random overlap” clouds are assigned an optical 
thickness based on an empirical relation between local temperature and optical thickness. 
The relation gives a thicker cloud when the temperature is higher and is expressed: 


0 

for 

T < 190.66 

(T~ 190.66) 2 *2.t10- 6 

for 

190.66 < T < 263.16 

6.95;rl0~ 3 * T — 1.82 

for 

263.16 < T < 273.38 

0.08 

for 

273.38 < T 


The relation for the range 190.66 < T < 263.16 is from Platt and Harshvardahn, 1988. The 
optical thickness for the longwave radiative feedback is assumed to be 75% of these values. 

The entire Moist Convective Processes Module is called with a frequency of 10 minutes. 
The cloud fraction values are time-averaged over the period between Radiation calls (every 
3 hours). Therefore, in a time- averaged sense, both random overlap and maximum overlap 
cloudiness can exist in a given grid-box. 


3.3 Turbulence 

The GEOS GCM turbulence parameterization consists of elements which handle vertical 
diffusion (Helfand and Labraga, 1988) and surface fluxes of heat, moisture and momentum 
(Helfand, et al, 1991, and Helfand and Schubert, 1994). The parameterization is invoked 
every 30 minutes, and employs a backward-implicit iterative time scheme with an internal 
time step of 5 minutes. The vertical regime is divided into a free atmosphere, a surface 
layer, and a viscous sub-layer above the surface roughness elements. The turbulent eddy 
fluxes are calculated using a variety of methods depending on the vertical location in the 
atmosphere. 

Turbulent eddy fluxes of momentum, heat and moisture in the surface layer are calculated 
using stability-dependant bulk formulae based on Monin-Obukhov similarity functions. For 
an unstable surface layer, the chosen stability functions are the KEYPS function (Panofsky, 
1973) for momentum, and its generalization for heat and moisture. The function for heat and 
moisture assures non- vanishing heat and moisture fluxes as the wind speed approaches zero. 
For a stable surface layer, the stability functions are those of Clarke (1970), slightly modified 
for the momemtum flux. The moisture flux also depends on a specified evapotranspiration 
coefficient, set to unity over oceans and dependant on the climatological ground wetness over 
land. The gradients in the viscous sublayer are based on Yaglom and Kader (1974). The 
surface roughness length over oceans is computed as an interpolation between the functions 
of Large and Pond (1981) for high winds and of Kondo (1975) for weak winds, and over 
land is specified from the climatology of Dorman and Sellers (1989). 
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Above the surface layer, turbulent fluxes of momentum., heat and moisture are calculated 
by the Level 2.5 Mellor-Yamada type closure scheme of Helfand and Labraga (1988), which 
predicts turbulent kinetic energy and determines the eddy transfer coefficients used for 
vertical diffusion. 


3.3.1 Planetary Boundary Layer 

The Planetary Boundary Layer (PBL) height is diagnosed by the parameterization as the 
level at which the turbulent kinetic energy is reduced to a tenth of its surface value, and 
generally occurs within the first two to four sigma levels above the surface (sigma equal 
0.994 to about sigma equal 0.875). 


4 Boundary Conditions 


The conditions at the lower boundary which are not explicitly predicted during execution 
must be specified from observed fields. Boundary condition data sets are available at present 
in monthly mean increments, at the model’s 4° x 5° and 2° x 2.5° resolutions, for either 
climatological or yearly varying conditions (beginning from 1979). Monthly mean values 
are interpolated during each model timestep to the current time. Future model versions 
will incorporate boundary conditions at higher spatial (1° x 1 ) and temporal resolutions. 


4.1 Orography 

Surface geopotential heights are provided from an averaging of the Navy 10 minute by 10 
minute dataset supplied by the National Center for Atmospheric Research (NCAR) to the 
model’s grid resolution. For the 2° x 2.5° GE03-1 GCM, the averaged topography was then 
passed through a Lanczos (1966) filter in both dimensions which removes the smallest scales 
while inhibiting Gibbs phenomena. Negative values in the topography resulting from the 
filtering procedure were not filled. 

In one dimension, we may define a cyclic function in x as: 

f(x) = ^ ^ {a k cos (kx) + b k sin (kx)) (21) 

1 k-1 

where N — and IM is the total number of points in the x direction. Defining Ax im’ 
we may define the average of f(x) over a 2 Ax region as: 
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Using equation (21) in equation (22) and integrating, we may write: 


a ° , 1 

/W = T + 2A^ 


N 

E 

k = 1 


OJb 


sin(A:a: / ) 

jfe 


/ 07+Aa: 
x—Ax 


b k 


cos(kx') 


ar-fAa;' 

x—Ax 


(23) 


or 


/(a?) = v + X Sm f ( a * cos (^) + sin(foc)) (24) 

2 *=i 

Thus, the Fourier wave amplitudes are simply modified by the Lanczos filter response func- 
tion s> - K |^y Y ^ ♦ This may be compared with an rath - order Shapiro (1970) filter response 
function, defined as 1 — sin m (^ £ ), shown in Figure 3. 


4.2 Land-Sea-Ice Surface Type 

The determination of the land or sea category of surface type was made from N CAR’s. 
10 minute by 10 minute Navy topography dataset, which includes information about the 
percentage of water-cover at any point. The data were averaged to the model’s 4° x 5° and 
2° x 2.5° grid resolutions. Any grid-box whose averaged water percentage was > 60% was 
defined as a water point. Permanent ice or glacier points for the 2° x 2.5° grid were identified 
from the climatology of Parkinson’s satellite-derived fields (NASA SP 459: Arctic Sea Ice, 
1973-1976, and NASA SP 489: Antartic Sea Ice, 1973-1976). The Land/Water masks and 
topography used for the 4° x 5° and 2° x 2.5°GEOS-l GCM grids are shown in Figure 4 
and Figure 5, respectively. 


4.3 Sea Surface Temperature 

Climatological sea surface temperatures for the 4° x 5° and 2° x 2.5° grid resolutions are 
interpolated from the data of Reynolds, 1988. Yearly varying monthly mean data for the 4° 
x 5° grid are from the AMIP specified monthly mean fields (Gates, 1992), and for the 2° x 
2.5° grid are interpolated from Reynolds’ (1988) data. Monthly mean data are interpolated 
linearly to the current time at every GCM timestep. 



Filter Response Functions 



Figure 3: Comparison between the Lanczos and mth - order Shapiro filter response functions 
for m = 2, 4, and 8. 
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Figure 4: GEOS-1 GCM Land/Water Mask and Topography used at 4° x 5° resolution. 
Light grey shading denotes land surfaces, while dark grey denotes permanent glacier. 
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Figure 5: GEOS-1 GCM Land/Water Mask and Topography used at 2° x 2.5° resolution. 
Light grey shading denotes land surfaces, while dark grey denotes permanent glacier. 
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4.4 Sea Surface Roughness 


Ocean surface roughness is computed iteratively as a function of wind stress in the Turbu- 
lence Scheme. 


4.5 Albedo 

Monthly mean albedos over land are data of Posey and Clapp (1964), with modifications as 
described in Kitzmiller (1979). Over oceans the albedo is specified as 0.07, and over ice (sea 
ice or glaciers) 0.80. Monthly mean values are linearly interpolated to the current time. 


4.6 Sea Ice 

Monthly mean sea ice extent were interpolated to the model’s 4° x 5° and 2° x 2.5° grid res- 
olutions from the AMIP monthly mean fields. Monthly mean data are linearly interpolated 
to the current time. Any sea ice is assumed to be three meters thick, and the conduction 
through sea ice is accounted for in the Turbulence Parameterization as part of the surface 
energy budget. 


4.7 Snow Cover 

Snow cover is assumed to occur when the ground temperature over the land is less than 
freezing in conjunction with climatological albedo values greater than 0.4. Surface conduc- 
tive characteristics are adjusted accordingly. 


4.8 Land Surface 

Monthly varying climatological roughness lengths are specified for each land surface vege- 
tation type from Dorman and Sellers (1989). A gridded surface vegetation type designation 
is from Sud, et al. (1990). Monthly mean values are interpolated linearly to the model’s 
current time. The soil wetness climatology and yearly varying fields are from the monthly 
estimates of Schemm, et al. (1992) based on the procedure developed by Mintz and Serafini 
(1984) using an inverted single layer “bucket” model together with observed fields of surface 
air temperature and precipitation. Values are interpolated linearly to the model’s current 
time. Ground temperature over land is predicted from a surface energy balance equation. 
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Part II 


GEOS-1 GCM: Diagnostics 

5 Introduction 


This section of the documentation describes the Diagnostics Utilities available within the 
Goddard Earth Observing Systems (GEOS) General Circulation Model (GCM). In addition 
to a description on how to set and extract diagnostic quantities, this document also provides 
a comprehensive list of all available diagnostic quantities and a short description of how 
they are computed. It should be noted that this document is not intended to be a complete 
documentation of the various physical parameterization schemes used in the GEOS GCM, 
and the reader should refer to original publications for further insight. 


6 Diagnostic Overview 


A large selection of model diagnostics is available in the GEOS GCM. At the time of this 
writing there are 119 different diagnostic quantities which can be enabled for an experi- 
ment. As a matter of philosophy, no diagnostic is enabled as default, thus each user must 
specify the exact diagnostic information required for an experiment. This is accomplished 
by enabling the specific diagnostic of interest cataloged in the Diagnostic Menu (see Section 
7). The Diagnostic Menu is a hard- wired enumeration of diagnostic quantities available 
within the GEOS GCM. Diagnostics are internally referred to by their associated number 
in the Diagnostic Menu. Once a diagnostic is enabled, the GEOS GCM will continually 
increment an array specifically allocated for that diagnostic whenever the associated process 
for the diagnostic is computed. Separate arrays are used both for the diagnostic quantity 
and its diagnostic counter which records how many times each diagnostic quantity has been 
computed. 

There are several utilities within the GEOS GCM available to users to enable, disable, 
clear, and retrieve model diagnostics, and may be called from any user-supplied application 
and/or output routine. The available utilities and the CALL sequences are listed below. 

SETDIAG: This subroutine enables a diagnostic from the Diagnostic Menu, meaning that 
space is allocated for the diagnostic and the model routines will increment the diagnostic 
value during execution. This routine is useful when called from either user application 
routines or user output routines, and is the underlying interface between the user and the 
desired diagnostic. The diagnostic is referenced by its diagnostic number from the menu, 
and its calling sequence is given by: 
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where 


CALL SETDIAG (NUM) 

NUM - Diagnostic number from menu 


GETDIAG: This subroutine retrieves the value of a model diagnostic. This routine is 
particulary useful when called from a user output routine, although it can be called from an 
application routine as well. This routine returns the time-averaged value of the diagnostic 
by dividing the current accumulated diagnostic value by its corresponding counter. This 
routine does not change the value of the diagnostic itself, that is, it does not replace the 
diagnostic with its time-average. The calling sequence for this routine is givin by: 


where 


CALL GETDIAG (LEV,NUM,QTMP,UNDEF) 

LEV = Sigma Level at which the diagnostic is desired 

NUM = Diagnostic number from menu 

QTMP = Time- Averaged Diagnostic Output 

UNDEF = Fill value to be used when diagnostic is undefined 


CLRDIAG: This subroutine initializes the values of model diagnostics to zero, and is 
particularly useful when called from user output routines to re- initialize diagnostics during 
the run. The calling sequence is: 


CALL CLRDIAG (NUM) 

where NUM = Diagnostic number from menu 


ZAPDIAG: This entry into subroutine SETDIAG disables model diagnostics, meaning 
that the diagnostic is no longer available to the user. The memory previously allocated to 
the diagnostic is released when ZAPDIAG is invoked. The calling sequence is given by. 

CALL ZAPDIAG (NUM) 

where NUM = Diagnostic number from menu 


DIAGSIZE: We end this section with a discussion on the manner in which computer 
memory is allocated for GEOS diagnostics. All GEOS GCM diagnostic quantities are 


24 



stored in the single diagnostic array QDIAG which is located in the DIAG COMMON, 
having the form: 


COMMON /DIAG/ NDIAG_MAX,QDIAG(IM,JM,1) 


where NDIAG-MAX is an Integer variable which should be set equal to the number of 
enabled diagnostics, and QDIAG is a three-dimensional array. The first two-dimensions 
of QDIAG correspond to the horizontal dimension of a given diagnostic, while the third 
dimension of QDIAG is used to identify specific diagnostic types. In order to minimize the 
maximum memory requirement used by the model, the GEOS GCM libraries are compiled 
with room for only one horizontal diagnostic array, as shown in the above example. In 
order for the User to enable more than 1 two-dimensional diagnostic, the size of the DIAG 
COMMON must be expanded to accomodate the desired diagnostics. This can be accom- 
plished by including the COMMON DIAG in the User’s Application program, with the 
QDIAG array properly dimensioned. Another method, currently used by DAO production 
routines, involves the use of the DIAGSIZE subroutine. Using the DIAGSIZE subroutine 
is a convenient method to allocate just enough permanent memory during a run to be 
used for diagnostic storage. All production DAO applications contain a call to subroutine 
DIAGSIZE, which expands the size of the DIAG COMMON to include the total number 
of surface and upper- air diagnostics the User has enabled. An example of the DIAGSIZE 
subroutine is shown below: 


subroutine diagsize 
parameter ( IM = 144 ) 
parameter ( JM = 91 ) 

parameter ( NLAY =20 ) 

parameter ( ndiags =11 ) 

parameter ( ndiagu = 4 ) 

common /diag/ ndiag_max,qdiag(im, jm,ndiags+ndiagu*nlay) 

ndiag_max = ndiags + ndiagu*nlay 

return 

end 
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7 GEOS GCM Diagnostic Menu 


N NAME UNITS LEVELS DESCRIPTION 


1 

UFLUX 

Newton /m 2 

1 

2 

VFLUX 

Newton jm 2 

1 

3 

HFLUX 

Watts/m 2 

1 

4 

EFLUX 

Watts jm 2 

1 

5 

QICE 

Watts jm 2 

1 

6 

RADLWG 

Watts jm 2 

1 

7 

RADSWG 

Watts/m 2 

1 

8 

RI 

deg /day 

NLAY 

9 

CT 

mj see 

1 

10 

CU 

mj sec 

1 

11 

ET 

m 2 / sec 

NLAY 

12 

EU 

m 2 j sec 

NLAY 

13 

TURBU 

mj sec /day 

NLAY 

14 

TURBV 

mj sec/ day 

NLAY 

15 

TURBT 

deg /day 

NLAY 

16 

TURBQ 

g/kg/day 

NLAY 

17 

MOISTT 

deg /day 

NLAY 

18 

MOISTQ 

g/kg/day 

NLAY 

19 

RADLW 

deg / day 

NLAY 

20 

RADSW 

deg j day 

NLAY 

21 

PREACC 

mm/ day 

1 

22 

PRECON 

mm/ day 

1 

23 

TUFLUX 

Newton jm 2 

NLAY 

24 

TVFLUX 

Newton /m 2 

NLAY 

25 

TTFLUX 

Watts/m 2 

NLAY 

26 

TQFLUX 

Watts jm 2 

NLAY 

27 

CN 

m/ sec 

1 

28 

WINDS 

mj sec 

1 

29 

DTSRF 

deg 

1 

30 

TG 

deg 

1 

31 

TS 

deg 

1 

32 

DTG 

deg 

1 


Surface U-Wind Stress on the atmosphere 
Surface V-Wind Stress on the atmosphere 
Surface Flux of Sensible Heat 
Surface Flux of Latent Heat 
Heat Conduction through Sea-Ice 
Net upward LW flux at the ground 
Net downward SW flux at the ground 
Richardson Number 
Surface Drag coefficient for T and Q 
Surface Drag coefficient for U and V 
Diffusivity coefficient for T and Q 
Diffusivity coefficient for U and V 
U-Momentum Changes due to Turbulence 
V-Momentum Changes due to Turbulence 
Temperature Changes due to Turbulence 
Specific Humidity Changes due to Turbulence 
Temperature Changes due to Moist Processes 
Specific Humidity Changes due to Moist Pro- 
cesses 

Net Longwave heating rate for each level 
Net Shortwave heating rate for each level 
Total Precipitation 
Convective Precipitation 
Turbulent Flux of U-Momentum 
Turbulent Flux of V-Momentum 
Turbulent Flux of Sensible Heat 
Turbulent Flux of Latent Heat 
Neutral Drag Coefficient 
Surface Wind Speed 

Air/Surface virtual temperature difference 
Ground temperature 

Surface air temperature (Adiabatic from low- 
est model layer) 

Ground temperature adjustment 
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N NAME UNITS LEVELS DESCRIPTION 


33 

QG 

g/kg 

1 

34 

QS 

g/kg 

1 

35 

TGRLW 

deg 

1 

36 

ST4 

Watts /m 2 

1 

37 

OLR ■ 

Watts /m 2 

1 

38 

OLRCLR 

W atts / m 2 

1 

39 

LWGCLR 

Watts jm 2 

1 

40 

LWCLR 

deg 1 day 

NLAY 

41 

TLW 

deg 

NLAY 

42 

SHLW 

gig 

NLAY 

43 

OZLW 

g/g 

NLAY 

44 

CLMOLW 

0-1 

NLAY 

45 

CLROLW 

0 - 1 

NLAY 

46 

CLMOSW 

0 - 1 

NLAY 

47 

CLROSW 

0 - 1 

NLAY 

48 

RADSWT 

W atts /m 2 

1 

49 




50 




51 




52 

EVAP 

mm /day 

1 

53 

DPDT 

mm /day 

1 

54 

ANALU 

m/ sec/ sec 

NLAY 

55 

ANALV 

mj sec/ sec 

NLAY 

56 

ANALT 

mb deg / mb K / sec 

NLAY 

57 

ANALQ 

mb * g / g j sec 

NLAY 

58 

OMEGA 

mb /day 

NLAY 


Ground specific humidity 
Saturation surface specific humidity 

Instantaneous ground temperature used as in- 
put to the Longwave radiation subroutine 
Upward Longwave flux at the ground (oT 4 ) 
Net upward Longwave flux at the top of the 
model 

Net upward clearsky Longwave flux at the top 
of the model 

Net upward clearsky Longwave flux at the 
ground 

Net clearsky Longwave heating rate for each 
level 

Instantaneous temperature used as input to 

the Longwave radiation subroutine 

Instantaneous specific humidity used as input 

to the Longwave radiation subroutine 

Instantaneous ozone used as input to the 

Longwave radiation subroutine 

Maximum overlap cloud fraction used in the 

Longwave radiation subroutine 

Random overlap cloud fraction used in the 

Longwave radiation subroutine 

Maximum overlap cloud fraction used in the 

Shortwave radiation subroutine 

Random overlap cloud fraction used in the 

Shortwave radiation subroutine 

Incident Shortwave radiation at the top of the 

atmosphere 

Disabled 

Disabled 

Disabled 

Surface evaporation 
Total surface pressure tendency 
Analysis u — Wind increment 
Analysis v — Wind increment 
Analysis n6 increment, 6 — ~ 

Analysis irq increment 
Vertical (Omega) velocity 
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N NAME UNITS LEVELS DESCRIPTION 


59 

DUDT 

m/ sec/day 

NLAY 

60 

DVDT 

m / sec f day 

NLAY 

61 

DTDT 

deg / day 

NLAY 

62 

DQDT 

g/kg/day 

NLAY 

63 

VORT 

sec -1 

NLAY 

64 




65 




66 




67 

USTAR 

m/ sec 

1 

68 

ZO 

m 

1 

69 

FRQTRB 

0-1 

NLAY-1 

70 

PBL 

mb 

1 

71 

SWCLR 

deg /day 

NLAY 

72 

OSR 

Watts /m 2 

1 

73 

OSRCLR 

Watts/m 2 

1 

74 

CLDMAS 

kg m/ sec 

NLAY 

75 

UAVE 

m/ sec 

NLAY 

76 

VAVE 

m/ sec 

NLAY 

77 

TAVE 

deg 

NLAY 

78 

QAVE 

(j/9 

NLAY 

79 




80 

PAVE 

mb 

1 

81 

QQAVE 

[mj sec) 2 

NLAY 

82 

SWGCLR 

Watts jm 2 

1 

83 

ANALP 

mb j sec 

1 

84 

SDIAG1 


1 

85 

SDIAG2 


1 

86 

UDIAG1 


NLAY 

87 

UDIAG2 


NLAY 

88 

DIABU 

m/ sec /day 

NLAY 

89 

DIABV 

m/ sec/day 

NLAY 

90 

DIABT 

deg /day 

NLAY 

91 

DIABQ 

g/kg/day 

NLAY 


Total U-Wind tendency 

Total V-Wind tendency 

Total Temperature tendency 

Total Specific Humidity tendency 

Relative Vorticity 

Disabled 

Disabled 

Disabled 

Surface USTAR wind 

Surface roughness 

Frequency of Turbulence 

Planetary Boundary Layer depth 

Net clear sky Shortwave heating rate for each 

level 

Net downward Shortwave flux at the top of 
the model 

Net downward clearsky Shortwave flux at the 
top of the model 
Convective cloud mass flux 
Time-averaged u — Wind 
Time-averaged v — Wind 
Time-averaged Temperature 
Time-averaged Specific Humidity 
Disabled 

Time-averaged p sur j - p top 
Time- averaged Turbulent KineticEnergy 
Net downward clearsky Shortwave flux at the 
ground 

Analysis Sur facePressure increment 
User- Defined Surface Diagnostic- 1 
User-Defined Surface Diagnostic-2 
User-Defined Upper- Air Diagnostic- 1 
User-Defined Upper-Air Diagnostic-2 
Total Diabatic forcing on u — Wind 
Total Diabatic forcing oni?- Wind 
Total Diabatic forcing on Temperature 
Total Diabatic forcing on Sped fic Humidity 
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N NAME 


UNITS 


LEVELS DESCRIPTION 


92 

93 

94 

95 

96 

97 

98 

99 

100 
101 


102 

VINTUQ 

m/ sec • g/kg 

1 

103 

VINTVQ 

m/ sec • g/kg 

1 

104 

VINTUT 

mj sec • deg 

1 

105 

VINTVT 

mj sec • deg 

1 

106 

CLDFRC 

0- 1 

1 

107 

QINT 

gm/cm 2 

1 

108 

U2M 

mj sec 

1 

109 

V2M 

mj sec 

1 

110 

T2M 

deg 

1 

111 

Q2M 

g/kg 

1 

112 

U10M 

m/sec 

1 

113 

V10M 

mj sec 

1 

114 

T10M 

deg 

1 

115 

Q10M 

g/kg 

1 

116 

DTRAIN 

kg * mj sec 

NLAY 

117 

QFILL 

g/kg/day 

NLAY 

118 

VINTQANA 

mm /day 

1 

119 

VINTQFIL 

mm/day 

1 


Disabled 

Disabled 

Disabled 

Disabled 

Disabled 

Disabled 

Disabled 

Disabled 

Disabled 

Disabled 

Vertically integrated u q 
Vertically integrated v q 
Vertically integrated uT 
Vertically integrated v T 
Total Cloud Fraction 
Precipitable water 
U-Wind at 2 meters 
V-Wind at 2 meters 
Temperature at 2 meters 
Specific Humidity at 2 meters 
U-Wind at 10 meters 
V-Wind at 10 meters 
Temperature at 10 meters 
Specific Humidity at 10 meters 
Convective Cloud Detrainment 
Filling of negative specific humidity 
Vertically integrated analysis specific humid- 
ity increment 

Vertically integrated QFILL 
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8 Diagnostic Description 


In this section we list and describe the diagnostic quantities available within the GEOS 
GCM. The diagnostics are listed in the order that they appear in the Diagnostic Menu, 
Section 7. In all cases, each diagnostic as currently archived on DAO production output 
datasets is time-averaged over its diagnostic output frequency: 

t=TTOT 

DIAGNOSTIC = E 

where TTOT — , NQDIAG is the output frequency of the diagnositc, and At is 

the timestep over which the diagnostic is updated. For further information on how to set 
the diagnostic output frequency NQDIAG, please see Part III, A User’s Guide. Also, note 
that several Diagnostic Index locations are currently disabled due to changes which have 
occured from previous GCM versions. 

A word of caution: The Arakawa C-grid stencil used in the GEOS GCM employs a 

“stretched” coordinate near the poles due to conservation requirements. As a result, a 
well-defined vorticity is is located at the poles, while the first mass point is located one 
A <f> away from the poles. Since all Physical (diabatic) processes of the model are computed 
at mass points, all Physics diagnostics are only truely defined at latitudes — 90 + A (j) to 
90 - A <j>. While an attempt was made to “diagnose” some pole value diagnostics through 
interpolation over the poles, others have been left unfilled. 


1) UFLUX Surface Zonal Wind Stress on the Atmosphere (Newton /m 2 ) 

The zonal wind stress is the turbulent flux of zonal momentum from the surface. See section 
3.3 for a description of the surface layer parameterization. 

UFLUX = -pC D W s u where : C D = C\ 

where p — the atmospheric density at the surface, Co is the surface drag coefficient, C u is the 
dimensionless surface exchange coefficient for momentum (see diagnostic number 10), W s 
is the magnitude of the surface layer wind, and u is the zonal wind in the lowest model layer. 


2) VFLUX Surface Meridional Wind Stress on the Atmosphere (Newton/m 2 ) 

The meridional wind stress is the turbulent flux of meridional momentum from the surface. 
See section 3.3 for a description of the surface layer parameterization. 

VFLUX = - pC D W s v where : C D = C 2 U 
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where p = the atmospheric density at the surface, C'd is the surface drag coefficient, C u is 
the dimensionless surface exchange coefficient for momentum (see diagnostic number 10), 
W s is the magnitude of the surface layer wind, and v is the meridional wind in the lowest 
model layer. 


3) HFLUX Surface Flux of Sensible Heat (Watts/m 2 ) 

The turbulent flux of sensible heat from the surface to the atmosphere is a function of the 
gradient of virtual potential temperature and the eddy exchange coefficient: 

HFLUX - P K pc p C H W s (0 sur f ace - 0nlay) where : Cjj = C u C t 

where p = the atmospheric density at the surface, c p is the specific heat of air, Ch is the 
dimensionless surface heat transfer coefficient, W s is the magnitude of the surface layer 
wind, C u is the dimensionless surface exchange coefficient for momentum (see diagnostic 
number 10), G't is the dimensionless surface exchange coefficient for heat and moisture (see 
diagnostic number 9), and 0 is the potential temperature at the surface and at the bottom 
(j-level. 


4) EFLUX Surface Flux of Latent Heat ( Watts/m 2 ) 

The turbulent flux of latent heat from the surface to the atmosphere is a function of the 
gradient of moisture, the potential evapotranspiration fraction and the eddy exchange co- 
efficient: 


EFLUX — pftLG jj\ / V s ((j[ sur f ace (Inlay ) where : Ch — C u Ct 

where p = the atmospheric density at the surface, (5 is the fraction of the potential evapo- 
transpiration actually evaporated, L is the latent heat of evaporation, Ch is the dimension- 
less surface heat transfer coefficient, W s is the magnitude of the surface layer wind, C u is 
the dimensionless surface exchange coefficient for momentum (see diagnostic number 10), 
Ct is the dimensionless surface exchange coefficient for heat and moisture (see diagnostic 
number 9), and q SU rface and (jnlay are the specific humidity at the surface and at the 
bottom (j-level, respectively. 


5) QICE Heat Conduction Through Sea Ice (Watts /m 2 ) 

Over sea ice there is an additional source of energy at the surface due to the heat conduction 
from the relatively warm ocean through the sea ice. The heat conduction through sea ice 
represents an additional energy source term for the ground temperature equation. 


QICE = - T g ) 
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where Cu is the thermal conductivity of ice, Hi is the ice thickness, assumed to be 3 m 
where sea ice is present, T{ is 273 degrees Kelvin, and T g is the temperature of the sea ice. 

NOTE: QICE is not available through model version 5.3, but is available in subsequent 
versions. 


6) RADLWG Net upward Longwave Flux at the surface (Watts /m 2 ) 


RADLWG 


rpNet 

*LW,nlay + 1 

^ •j' j | 

^LW ) nlay -\- 1 LW, nlay+1 


where nlay+1 indicates the lowest model edge-level, or p — p SU rf • F'lw * s u P war d Long- 
wave flux and F^ w is the downward Longwave flux. 


7) RADSWG Net downard shortwave Flux at the surface (Watts/m 2 ) 


RADSWG 


rpN et 

^SW, nlay+1 

-t T 

^SW, nlay+1 ^SW, nlay+1 


where nlay+1 indicates the lowest model edge-level, or p = p SU rf- Fg W downward 

Shortwave flux and F$ w is the upward Shortwave flux. 


8) RI Richardson Number ( dimensionless ) 

The non-dimensional stability indicator is the ratio of the buoyancy to the shear: 


RI - 


a dd v 
8 V dz 


( i ) 2 + (§?) 2 


_ de„ dP K 
dz dz 

df) 2 +(i ) 2 


where we used the hydrostatic equation: 


OP* 


— CpO v 


Negative values indicate unstable buoyancy AND shear, small positive values (< 0.4) 
indicate dominantly unstable shear, and large positive values indicate dominantly stable 
stratification. 
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9) CT Surface Exchange Coefficient for Temperature and Moisture ( dimensionless ) 


The surface exchange coefficient is obtained from the similarity functions for the stability 
dependant flux profile relationships: 


CT = = £ w ' q ‘ = k 

u*A0 u*Aq {i>h + i>g) 

where ^ is the surface layer non-dimensional temperature change and ip a is the viscous 
sublayer non-dimensional temperature or moisture change: 


*l>h ~ f Y dC> 
k o C 


and 


, 0.55(-Pr 2 / 3 - 0.2) t n1/2 

i >9 = ^172 «♦,«,) ' 


and: /iq = 30^o with a maximum value over land of 0.01 


(j>h is the similarity function of which expresses the stability dependance of the temper- 
ature and moisture gradients, specified differently for stable and unstable layers according 
to Helfand and Schubert, 1993. k is the Von Karman constant, ( is the non-dimensional 
stability parameter, Pr is the Prandtl number for air, v is the molecular viscosity, zq is the 
surface roughness length, u * is the surface stress velocity (see diagnostic number 67), and 
the subscript ref refers to a reference value. 


10) CU Surface Exchange Coefficient for Momentum ( dimensionless ) 


The surface exchange coefficient is obtained from the similarity functions for the stability 
dependant flux profile relationships: 


CU 


w* 

w s 


k 

iftrn 


where is the surface layer non-dimensional wind shear: 



d> m is the similarity function of £, which expresses the stability dependance of the temper- 
ature and moisture gradients, specified differently for stable and unstable layers according 
to Helfand and Schubert, 1993. k is the Von Karman constant, ( is the non-dimensional 
stability parameter, w* is the surface stress velocity (see diagnostic number 67), and W s is 
the magnitude of the surface layer wind. 


11) ET Diffusivity Coefficient for Temperature and Moisture (in 2 /sec) 
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In the level 2.5 version of the Mellor-Yamada (1974) hierarchy, the turbulent heat or mois- 
ture flux for the atmosphere above the surface layer can be expressed as a turbulent diffusion 
coefficient Kh times the negative of the gradient of potential temperature or moisture. In 
the Helfand and Labraga (1988) adaptation of this closure, Kh takes the form: 


ET - K h 



qtSn(GM,GH) for decaying turbulence 
^ i SH(GM e ,GH e ) for growing turbulence 


where q is the turbulent velocity, or y/2 * turbulent kinetic energy, q e is the turbulence 
velocity derived from the more simple level 2.0 model, which describes equilibrium turbu- 
lence, t is the master length scale related to the layer depth, Sh is a function of Gh and 
Gmi the dimensionless buoyancy and wind shear parameters, respectively, or a function of 
Gh c and (?M e , the equilibrium dimensionless buoyancy and wind shear parameters. Both 
Gh and Gm, and their equilibrium values Gn e and Gm c , are functions of the Richardson 
number. 


For the detailed equations and derivations of the modified level 2.5 closure scheme, see 
Helfand and Labraga, 1988. 

In the surface layer, ET is the exchange coefficient for heat and moisture, in units of m/sec , 
given by: 

ETnlaY - C t * = Cf{W 8 

where Ct is the dimensionless exchange coefficient for heat and moisture from the surface 
layer similarity functions (see diagnostic number 9), u * is the surface friction velocity (see 
diagnostic number 67), Ch is the heat transfer coefficient, and W s is the magnitude of the 
surface layer wind. 


12) EU Diffusivity Coefficient for Momentum ( m 2 /sec ) 


In the level 2.5 version of the Mellor-Yamada (1974) hierarchy, the turbulent heat momen- 
tum flux for the atmosphere above the surface layer can be expressed as a turbulent diffusion 
coefficient I( m times the negative of the gradient of the u-wind. In the Helfand and Labraga 
(1988) adaptation of this closure, K m takes the form: 


EU - K m = 


( u'uj ') 
du 

dz 


qt Sm{Gmi Gh) for decaying turbulence 
^ l Sm ( Gm £ , G H e ) for growing turbulence 


where q is the turbulent velocity, or \J2 * turbulent kinetic energy , q e is the turbulence 
velocity derived from the more simple level 2.0 model, which describes equilibrium turbu- 
lence, l is the master length scale related to the layer depth, Sm is a function of Gh and 
Gm-, the dimensionless buoyancy and wind shear parameters, respectively, or a function of 
Gh £ and Gm g -, the equilibrium dimensionless buoyancy and wind shear parameters. Both 
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Gh and Gm, and their equilibrium values Gn e and G'M e > are functions of the Richardson 
number. 

For the detailed equations and derivations of the modified level 2.5 closure scheme, see 
Helfand and Labraga, 1988. 

In the surface layer, EU is the exchange coefficient for momentum, in units of m/sec, given 
by: 

EUnlay = C u * u* = CdW s 

where C u is the dimensionless exchange coefficient for momentum from the surface layer 
similarity functions (see diagnostic number 10), u * is the surface friction velocity (see di- 
agnostic number 67), Cd is the surface drag coefficient, and W s is the magnitude of the 
surface layer wind. 


13) TURBU Zonal U-Momentum changes due to Turbulence (m/ sec/ day) 


The tendency of U-Momentum due to turbulence is written: 

TT _._, TT du d —r-,\ & /T , ® u \ 

TURBU = — = = TrAmd 

Ot turb OZ OZ OZ 


The Helfand and Labraga level 2.5 scheme models the turbulent flux of u-momentum in 
terms of K m , and the equation has the form of a diffusion equation. 

14) TURBV Meridional V-Momentum changes due to Turbulence (m/ sec/ day) 


The tendency of V-Momentum due to turbulence is written: 


TURBV = 


dv 


d 


d 


dv 


o, = 7 r(-vV) = «-(Am7T“) 

Ot turb OZ OZ OZ 


The Helfand and Labraga level 2.5 scheme models the turbulent flux of v-momentum in 
terms of K m , and the equation has the form of a diffusion equation. 


15) TURBT Temperature changes due to Turbulence (deg /day) 
The tendency of temperature due to turbulence is written: 


The Helfand and Labraga level 2.5 scheme models the turbulent flux of temperature in 
terms of AT, and the equation has the form of a diffusion equation. 
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16) TURBQ Specific Humidity changes due to Turbulence (g/ kg /day) 


The tendency of specific humidity due to turbulence is written: 

TURBQ = ^ 

dtturb dz dz dz 


The Helfand and Labraga level 2.5 scheme models the turbulent flux of temperature in 
terms of Kh, and the equation has the form of a diffusion equation. 


17) MOISTT Temperature Changes Due to Moist Processes (deg /day) 


where: 


and 


The subscript c refers to convective processes, while the subscript Is refers to large scale 
precipitation processes, or supersaturation rain. The summation refers to contributions 
from each cloud type called by RAS. The dry static energy is given as s, the convective 
cloud base mass flux is given as m B , and the cloud entrainment is given as 77 , which are 
explicitly defined in Section 3.1.1, the description of the convective parameterization. The 
fractional adjustment, or relaxation parameter, for each cloud type is given as a, while R 
is the rain re-evaporation adjustment. 



18) MOISTQ Specific Humidity Changes Due to Moist Processes ( g/kg/day ) 


MOISTQ = 


dq 

dt 



Is 


where: 


and 


Oq 
dt c 




and 


dq 

dt 


= (9* - 9) 


r* = 99 


ds 

dp 


and 


I \ = 9V 


dh 

dp 


The subscript c refers to convective processes, while the subscript Is refers to large scale 
precipitation processes, or supersaturation rain. The summation refers to contributions 
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from each cloud type called by RAS. The dry static energy is given as s, the moist static 
energy is given as h, the convective cloud base mass flux is given as mg, and the cloud 
entrainment is given as t 7 , which are explicitly defined in Section 3.1.1, the description of 
the convective parameterization. The fractional adjustment, or relaxation parameter, for 
each cloud type is given as a, while R is the rain re-evaporation adjustment. 


19) RADLW Heating Rate due to Longwave Radiation (deg J day) 

The net longwave heating rate is calculated as the vertical divergence of the net terrestrial 
radiative fluxes. Both the clear-sky and cloudy-sky longwave fluxes are computed within the 
longwave routine. The subroutine calculates the clear-sky flux, F^ rsky , first. For a given 
cloud fraction, the clear line-of-sight probability C(p,p') is computed from the current level 
pressure p to the model top pressure, p 1 = p top , and the model surface pressure, p' = p S urf-> 
for the upward and downward radiative fluxes, (see Section 3.2.3). The cloudy-sky flux is 
then obtained as: 

Flw = C(p, /)■!&*•, 

Finally, the net longwave heating rate is calculated as the vertical divergence of the net 
terrestrial radiative fluxes: 

dpCpT 0 j^net 

~dT~~Vz FLW ’ 

or 

RADLW = 

Cp'K 0(7 

where g is the accelation due to gravity, c p is the heat capacity of air at constant pressure, 
and 

rpNET r>t rpi 

r LW ~ l LW “ l LW 


20) RADSW Heating Rate due to Shortwave Radiation (deg /day) 

The net Shortwave heating rate is calculated as the vertical divergence of the net solar 
radiative fluxes. The clear-sky and cloudy-sky shortwave fluxes are calculated separately. 
For the clear-sky case, the shortwave fluxes and heating rates are computed with both 
CLMO (maximum overlap cloud fraction) and CLRO (random overlap cloud fraction) set to 
zero (see Section 3.2.3). The shortwave routine is then called a second time, for the cloudy- 
sky case, with the true time-averaged cloud fractions CLMO and CLRO being used. In all 
cases, a normalized incident shortwave flux is used as input at the top of the atmosphere. 
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The heating rate due to Shortwave Radiation under cloudy skies is defined as: 


F (cloudy) fjp ■ RADSWT, 


RADSW = ~ - S- F(cloudy)sw T • RADSWT. 
c k 7 t oa 


where g is the accelation due to gravity, c p is the heat capacity of air at constant pres- 
sure, RADSWT is the true incident shortwave radiation at the top of the atmosphere (See 
Diagnostic #48), and 


F(cloudy)sw = F(doudy)\ w - F(cloudy) l sw 


21) PREACC Total (Large-scale + Convective) Accumulated Precipition (mm/ day) 


For a change in specific humidity due to moist processes, Aq mo i s t, the vertical integral or 
total precipitable amount is given by: 

f top , f top . dp it f 1 . , 

PREACC = / pAq mo i s tdz / Aq mo { s t / Aq mo i s td<7 

Jsurf Jsurf 9 9 J 0 


A precipitation rate is defined as the vertically integrated moisture adjustment per Moist 
Processes time step, scaled to mm/ day. 


22) PRECON Convective Precipition (mm/ day) 

For a change in specific humidity due to sub-grid scale cumulus convective processes, Aq cum , 
the vertical integral or total precipitable amount is given by: 

f to P , f t0 P A dp 7T f 1 . . 

PRECON = / pAq cum dz — — I Aq cum — ~~ Aq cum d(j 
Jsurf Jsurf 9 9 J 0 


A precipitation rate is defined as the vertically integrated moisture adjustment per Moist 
Processes time step, scaled to mm /day. 
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23) TUFLUX Turbulent Flux of U-Momentum (Newton /m 2 ) 


The turbulent flux of u-momentum is calculated for diagnostic purposes only from the 
eddy coefficient for momentum: 


A j r 

TUFLUX = = p (-K m —) 


where p is the air density, and K m is the eddy coefficient. 


24) TVFLUX Turbulent Flux of V-Momentum (Newton/ m 2 ) 

The turbulent flux of v-momentum is calculated for diagnostic purposes only from the 
eddy coefficient for momentum: 


TVFLUX 


P(V'W') = p(-Km^) 


where p is the air density, and K m is the eddy coefficient. 


25) TTFLUX Turbulent Flux of Sensible Heat (Watts/m 2 ) 

The turbulent flux of sensible heat is calculated for diagnostic purposes only from the 
eddy coefficient for heat and moisture: 


TTFLUX = CppP K (w , 0') = c pP P K (-K h -£) 


where p is the air density, and Kh is the eddy coefficient. 


26) TQFLUX Turbulent Flux of Latent Heat (Watts/m 2 ) 

The turbulent flux of latent heat is calculated for diagnostic purposes only from the eddy 
coefficient for heat and moisture: 

TQFLUX = Lp(W) = Lp(-Kh^) 


where p is the air density, and Kh is the eddy coefficient. 
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27) CN Neutral Drag Coefficient ( dimensionless ) 

The drag coefficient for momentum obtained by assuming a neutrally stable surface layer: 



where k is the Von Karman constant, h is the height of the surface layer, and z 0 is the 
surface roughness. 

NOTE: CN is not available through model version 5.3, but is available in subsequent ver- 
sions. 

28) WINDS Surface Wind Speed (meter /sec) 

The surface wind speed is calculated for the last internal turbulence time step. 

WINDS = , Ju% LA y + V% LA y 

where the subscript NLAY refers to the lowest model a - level. 

29) DTSRF Air/Surface Virtual Temperature Difference ( deg K) 

The air/surface virtual temperature difference measures the stability of the surface layer: 

DTSRF = (OvNLAY+1 - Ovnlay)PsutJ 

where 

OvNLAY+1 “ — ( 1 4‘*609^ f /VLAV+l) and qNLAY+l = <!NLAY+P(<l*( T g> F s)-qNLAY) 

r surf 

(3 is the surface potential evapotranspiration coefficient (0 = 1 over oceans), q*(T g , P s ) is the 
saturation specific humidity at the ground temperature and surface pressure, level NLAY 
refers to the lowest model level and level N LAY + 1 refers to the surface. 

30) TG Ground Temperature ( deg Ii) 

The ground temperature equation is solved as part of the turbulence package using a back- 
ward implicit time differencing scheme: 

dT 

TG is obtained from : C g —~ — R S w — Riw + Qice — H — LE 
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where R sw is the net surface downward shortwave radiative flux, Ri w is the net surface 
upward longwave radiative flux, Qi ce is the heat conduction through sea ice, H is the upward 
sensible heat flux, LE is the upward latent heat flux, and C g is the total heat capacity of 
the ground. C g is obtained by solving a heat diffusion equation for the penetration of the 
diurnal cycle into the ground (Blackadar, 1977), and is given by: 

^ = \/(0.386 + 0.536 W + 0.15W r2 )2a:10- 3 ^^: . 

2o? V 27 t 

Here, the thermal conductivity, A, is equal to 2x 10 -3 the angular velocity of the 

earth, u>, is written as 86400 sec/ day divided by 2n radians/ day, and the expression for 
C s , the heat capacity per unit volume at the surface, is a function of the ground wetness, W. 



31) TS Surface Temperature ( deg K) 

The surface temperature estimate is made by assuming that the model’s lowest layer is 
well-mixed, and therefore that 9 is constant in that layer. The surface temperature is 
therefore: 

TS = d^LAYPsurf 


32) DTG Surface Temperature Adjustment ( deg K) 

The change in surface temperature from one turbulence time step to the next, solved using 
the Ground Temperature Equation (see diagnostic number 30) is calculated: 

DTG = T g n - T g n ~ l 

where superscript n refers to the new, updated time level, and the superscript n — 1 refers 
to the value at the previous turbulence time level. 


33) QG Ground Specific Humidity ( g/kg ) 

The ground specific humidity is obtained by interpolating between the specific humidity at 
the lowest model level and the specific humidity of a saturated ground. The interpolation 
is performed using the potential evapotranspiration function: 

QG = qNLAY+l = QNLAY + P(q*(T g ,P s ) - qNLAY ) 


where j3 is the surface potential evapotranspiration coefficient ((3=1 over oceans), and 
q*(T g ,P s ) is the saturation specific humidity at the ground temperature and surface pres- 
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sure. 


34) QS Saturation Surface Specific Humidity (g/kg) 

The surface saturation specific humidity is the saturation specific humidity at the ground 
temprature and surface pressure: 

QS = q*{T gi P.) 

35) TGRLW Instantaneous ground temperature used as input to the Longwave 
radiation subroutine (deg) 

TGRLW = T g {\ y (j),n) 

where T g is the model ground temperature at the current time step n. 

36) ST4 Upward Longwave flux at the surface (Watts/m 2 ) 

ST4 = <tT 4 

where a is the Stefan- Boltzmann constant and T is the temperature. 

37) OLR Net upward Longwave flux at p — p top (Watts/ m 2 ) 

OLR = F^lv 

where top indicates the top of the first model layer. In the GEOS-1 GGM, pt op = 10 mb. 

38) OLRCLR Net upward clearsky Longwave flux at p- p top (Watts /m 2 ) 

OLRCLR = F(clearsky)%% T top 

where top indicates the top of the first model layer. In the GEOS-1 GGM, p top — 10 mb. 

39) LWGCLR Net upward clearsky Longwave flux at the surface (Watts /m 2 ) 

LWGCLR = F(clearsky)Lw,niay+i 

= F(clearsky)[ Wnlay+l - F(clearsky) l LWnlay+l 
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f 

where nlay-fl indicates the lowest model edge-level, or p = p SU rf • F(clearsky)' LW is the 
upward clearsky Longwave flux and the F(clearsky)\ jW is the downward clearsky Longwave 
flux. 


40) LWCLR Heating Rate due to Clearsky Longwave Radiation (deg /day) 

The net longwave heating rate is calculated as the vertical divergence of the net terrestrial 
radiative fluxes. Both the clear-sky and cloudy-sky longwave fluxes are computed within the 
longwave routine. The subroutine calculates the clear-sky flux, F^^ rsky , first. For a given 
cloud fraction, the clear line-of-sight probability C(p,p') is computed from the current level 
pressure p to the model top pressure, p’ — p top , and the model surface pressure, p' — p SU rf , 
for the upward and downward radiative fluxes, (see Section 3.2.3). The cloudy-sky flux is 
then obtained as: 


Flw = C(p,p') ■ 


Thus, LWCLR is defined as the net longwave heating rate due to the vertical divergence 
of the clear-sky longwave radiative flux: 


dpc p T 

dt clearsky 


-^-F(clearsky) 1 l § r , 


r\ 

LWCLR = F( clear sky )^ T . 

C p 7T OO 

where g is the accelation due to gravity, c p is the heat capacity of air at constant pressure, 
and 

F (clecirsky)Lw — F(clearsky)\ w — F(clearsky )^ LW 


41) TLW Instantaneous temperature used as input to the Longwave radiation 
subroutine (deg) 

TLW = T(\,4>,o,n) 

where T is the model temperature at the current time step n. 


42) SHLW Instantaneous specific humidity used as input to the Longwave ra- 
diation subroutine (kg/kg) 


SHLW = q(\,(f>,o,n) 
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where q is the model specific humidity at the current time step n . 


43) QZLW Instantaneous ozone used as input to the Longwave radiation sub- 
routine (kg/kg) 

OZLW = OZ(A,<£,<r,n) 

where OZ is the interpolated ozone data set from the climatological monthly mean zonally 
averaged ozone data set. 


44) CLMOLW Maximum Overlap cloud fraction used in LW Radiation (0-1) 

CLMOLW is the time- averaged maximum overlap cloud fraction that has been filled by the 
Relaxed Arakawa/ Schubert Convection scheme and will be used in the Longwave Radiation 
algorithm. These are convective clouds whose radiative characteristics are assumed to be 
correlated in the vertical. For a complete description of cloud/radiative interactions, see 
Section 3.2.3. 

CLMOLW = CLMO ras ,lw( A,<M) 

Note, in Version 1 of the GEOS GCM, shortwave and longwave cloud fields are computed 
identically, ie. 

CLMOLW = CLMOSW 


45) CLROLW Random Overlap cloud fraction used in LW Radiation (0 — 1) 

CLROLW is the time-averaged random overlap cloud fraction that has been filled by 
the Relaxed Arakawa/ Schubert and Large-scale Convection schemes and will be used in the 
Longwave Radiation algorithm. These are convective and large-scale clouds whose radiative 
characteristics are not assumed to be correlated in the vertical. For a complete description 
of cloud/radiative interactions, see Section 3.2.3. 

CLROLW = C LRO RAS,LargeScale,L,w{^-> a ) 


Note, in Version 1 of the GEOS GCM, shortwave and longwave cloud fields are computed 
identically, ie. 

CLROLW = CLROSW 


46) CLMOSW Maximum Overlap cloud fraction used in SW Radiation (0 — 1) 
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CLMOSW is the time-averaged maximum overlap cloud fraction that has been filled by the 
Relaxed Arakawa/ Schubert Convection scheme and will be used in the Shortwave Radiation 
algorithm. These are convective clouds whose radiative characteristics are assumed to be 
correlated in the vertical. For a complete description of cloud/radiative interactions, see 
Section 3.2.3. 

CLMOSW = CLMORASjwiKho) 

Note, in Version 1 of the GEOS GCM, shortwave and longwave cloud fields are computed 
identically, ie. 

CLMOSW = CLMOLW 


47) CLRQSW Random Overlap cloud fraction used in SW Radiation (0 ~ 1) 

CLROSW is the time-averaged random overlap cloud fraction that has been filled by 
the Relaxed Arakawa/ Schubert and Large-scale Convection schemes and will be used in the 
Shortwave Radiation algorithm. These are convective and large-scale clouds whose radiative 
characteristics are not assumed to be correlated in the vertical. For a complete description 
of cloud/radiative interactions, see Section 3.2.3. 

CLROSW = C LRORAS,LargeScale,S\v(^, Cr) 


Note, in Version 1 of the GEOS GCM, shortwave and longwave cloud fields are computed 
identically, ie. 


CLROSW = CLROLW 


48) RADSWT Incident Shortwave radiation at the top of the atmosphere (Watts /m 2 ) 

RADS WT = || • cos(f> z 

Rl 

where So, is the extra-terrestial solar contant, R a is the earth-sun distance in Astronomical 
Units, and cos<j> z is the cosine of the zenith angle. It should be noted that RADSWT, 
as well as OSR and OSRCLR, are calculated at the top of the atmosphere (p=0 mb). 
However, the OLR and OLRCLR diagnostics are currently calculated at p — p top (10 mb 
for the 20-level GEOS-1 GCM). 


49) Disabled 
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50) Disabled 


51) Disabled 


52) EVAP Surface Evaporation (mm /day) 

The surface evaporation is a function of the gradient of moisture, the potential evapotran- 
spiration fraction and the eddy exchange coefficient: 

EVAP = p(3Kh(q S urface ~ (INLAY ) 

where p = the atmospheric density at the surface, / 3 is the fraction of the potential evapo- 
transpiration actually evaporated ((3 = 1 over oceans), Kh is the turbulent eddy exchange 
coefficient for heat and moisture at the surface in m/sec and qsurface and qNLAY are 
the specific humidity at the surface (see diagnostic number 34) and at the bottom cr-level, 
respectively. 


53) DPDT Total Surface Pressure Tendency (mb j day) 


DPDT is the total time-tendency of the surface pressure due to Hydrodynamic and Analysis 
forcing. There are no surface pressure changes due to physical diabatic processes, such as 
changes in total precipitable water amounts. 


DPDT = 


Qtt dir 

dt Dynamics dt Analysis 


where 7 r = p sur f - p top - In forecast simulations ^ Analysis - 0. 


54) AN ALU Analysis Zonal U-Wind Tendency (m/ sec/ sec) 

ANALU is the Zonal U-Wind tendency created during the assimilation cycle which is used 
in the Incremental Analysis Updating (IAU) procedure. 

ANALU = (u A (\,<t>,a,t)-u FG (\,<t>,<7,tj)/At A nalysis 

where u A corresponds to the u-wind after the analysis, u FG corresponds to the u-wind first 
guess, and A t An alysis is the time-step between analyses (6 hours). Note, ANALU is stored 
on the staggered C-grid. 


55) ANALV Analysis Meridional V-Wind Tendency (m/ sec/ sec) 
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ANALV is the Meridional V-Wind tendency created during the assimilation cycle which 
is used in the Incremental Analysis Updating (IAU) procedure. 

ANALV = (v A (\,<j),(T,t)- V FG (\,(l>,Cr,tf) /At Analysis 

where v A corresponds to the v-wind after the analysis, v FG corresponds to the v-wind first 
guess, and At Analysis is the time-step between analyses (6 hours). Note, ANALV is stored 
on the staggered C-grid. 


56) AN ALT Analysis Temperature Tendency (mb deg / mb K / sec) 

AN ALT is the pressure weighted Potential Temperature tendency created during the as- 
similation cycle which is used in the Incremental Analysis Updating (IAU) procedure. 

ANALT = [(xO) A (\,<f),(T,t) - (ir0) FG (\,4>, a, t)) /At Analysis 

where (i rO) A corresponds to the pressure weighted Potential Temperature after the analysis, 
( 7 t6) fg corresponds to the pressure weighted Potential Temperature first guess, 0 — T/p K , 
and At Analysis is the time-step between analyses (6 hours). 


57) ANALQ Analysis Temperature Tendency (mb-g/g/sec) 

ANALQ is the pressure weighted Specific Humidity tendency created during the assimila- 
tion cycle which is used in the Incremental Analysis Updating (IAU) procedure. 

ANALQ = ((xq) A (\, 0, <7, t ) - ( irq) FG {\ , 4 >, (7, t)) / At Analysis 

where (7 rq) A corresponds to the pressure weighted Specific Humidity after the analysis, 
(7 rq) FG corresponds to the pressure weighted Specific Humidity first guess, and At Analysis 
is the time-step between analyses (6 hours). 


58) OMEGA Pressure Coordinate Vertical Velocity (mb/ day) 


OMEGA, which is the change in pressure in time following the parcel, is diagnostically 
obtained by integrating the GCM’s continuity equation. We have 


since p — ira + prop- Now, 


OMEGA = % = (7T<j + 7T<j) 
dt v ' 


7T 

dir 

m 


dir _ 

- + v .w 
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Integrating the continuity equation from the top of the model to a given model edge-level, 
< 7 e , and using the boundary condition b \ a =o= 0, we have 

r°e. dir 

irb e = - J ^ V • (rV)da - o e — 

Combining the above equations, we have 

OMEGA = - / 6 V • (7rV)d<7 + a e V • Vtt 
Jo 

Note that in this formulation, OMEGA is computed at GCM model edge-levels. OMEGA 
at the model top edge is taken to be 0, thus the first computed value of OMEGA is at 
sigma edge 2, or the bottom edge of the top model sigma layer. 


59) DUDT Total Zonal U-Wind Tendency (m/ sec/ day) 


DUDT is the total time-tendency of the Zonal U-Wind due to Hydrodynamic, Diabatic, 
and Analysis forcing. 


DUDT = 


du du da dii 

dt Dynamics dt Moist dt Turbulence dt Analysis 


Note, in Version 1 of the GEOS GCM there are no u-wind tendencies due to cumulus fric- 
tion,’(ie. §f Mojs( = 0). Also, in forecast simulations ^ Anatysis = 0. 


60) DVDT Total Zonal V-Wind Tendency (m/ sec i /day) 

DVDT is the total time-tendency of the Meridional V-Wind due to Hydrodynamic, Dia- 
batic, and Analysis forcing. 

dv dv +dv_ ^ dv 

dt Dynamics dt Moist dt Turbulence dt Analysis 

Note, in Version 1 of the GEOS GCM there are no v-wind tendencies due to cumulus fric- 
tion, (ie. % Moist = 0). Also, in forecast simulations %} Analysis = 0. 


61) DTDT Total Temperature Tendency (deg /day) 

DTDT is the total time-tendency of Temperature due to Hydrodynamic, Diabatic, and 


Analysis forcing. 

DTDT = 


+ 


dT + ^L + — 

dt Dynamics dt MoistProcesses dt ShortwaveRadiation 

dT dT , dT 

_j -| 

dt LongwaveRadiation dt Turbulence dt Analysis 
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Note, this diagnostic is obtained through differentiation by parts and re-arranging the model 
prognostic equations. Since the thermodynamic equation is written in the mass-weighted 
flux form, we have 


dr 

dt 


P_ 
7 r 


dirO ^dn 

+ 


dt 


G'K K 


dt L p 


- 1 


where 0 — T /p K . Also, in forecast simulations, ^ 


Analysis 


0. 


62) DQDT Total Specific Humidity Tendency (g/ kg /day) 

DQDT is the total time-tendency of Specific Humidity due to Hydrodynamic, Diabatic, 
and Analysis forcing. 


DQDT = 


dq +dq_ + dq + dq 

dt Dynamics dt MoistProcesses dt Turbulence dt Analysis 


Note, this diagnostic is obtained through differentiation by parts and re-arranging the model 
prognostic equations. Since the moisture equation is written in the mass- weighted flux form, 

WehaTO dq (d* q dir\ 

* = Ur-W* 


Also, in forecast simulations, 


- o. 

dt Analysis 


63) VORT Relative Vorticity (sec x ) 

The relative vorticity used in the formulation of the potential energy and enstrophy conserv- 
ing scheme described in Section 2 is stored as the diagnostic VORT. The relative vorticity 
is given as: 

1 T dv d(u cos <f>Y 

VORT = ttt ~ ~ o ; 

a cos cp LoA o<p J 

This diagnostic is stored on the Arakawa C-grid used in the GEOS GCM at the “vorticity- 
point” locations, ie., centered between the u-points and v-points, including values well 
defined at both poles. 

64) Disabled 

65) Disabled 
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66 ) Disabled 


67) USTAR Surface-Stress Velocity (m/sec) 

The surface stress velocity, or the friction velocity, is the wind speed at the surface layer 
top impeded by the surface drag: 

k 

USTAR = C U W S where : C u = — - 

m 

C u is the non-dimensional surface drag coefficient (see diagnostic number 10), and W s is 
the surface wind speed (see diagnostic number 28). 

68 ) ZO Surface Roughness Length (m) 

Over the land surface, the surface roughness length is interpolated to the local time from 
the monthly mean data of Dorman and Sellers (1989). Over the ocean, the roughness length 
is a function of the surface- stress velocity, ti*. 

ZO — Ci'n* + C 2 U 2 + C 3 W* + C 4 4 


where the constants are chosen to interpolate between the reciprocal relation of Kondo(1975) 
for weak winds, and the piecewise linear relation of Large and Pond(1981) for moderate to 
large winds. 


69) FRQTRB Frequency of Turbulence (0-1) 

The fraction of time when turbulence is present is defined as the fraction of time when the 
turbulent kinetic energy exceeds some minimum value, defined here to be 0.005 m 2 /sec 2 . 
When this criterion is met, a counter is incremented. The fraction over the averaging in- 
terval is reported. 


70) PBL Planetary Boundary Layer Depth (mb) 

The depth of the PBL is defined by the turbulence parameterization to be the depth at 
which the turbulent kinetic energy reduces to ten percent of its surface value. 


PBL = Ppbl Psurface 

where P PB l is the pressure in mb at which the turbulent kinetic energy reaches one tenth 
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of its surface value, and P s is the surface pressure. 


71) SWCLR Clear sky Heating Rate due to Shortwave Radiation (deg / day) 

The net Shortwave heating rate is calculated as the vertical divergence of the net solar 
radiative fluxes. The clear-sky and cloudy-sky shortwave fluxes are calculated separately. 
For the clear-sky case, the shortwave fluxes and heating rates are computed with both 
CLMO (maximum overlap cloud fraction) and CLRO (random overlap cloud fraction) set to 
zero (see Section 3.2.3). The shortwave routine is then called a second time, for the cloudy- 
sky case, with the true time-averaged cloud fractions CLMO and CLRO being used. In all 
cases, a normalized incident shortwave flux is used as input at the top of the atmosphere. 

The heating rate due to Shortwave Radiation under clear skies is defined as: 

^ = -j-F(dear)^ • RADSWT, 
or 

A 

SWCLR = —^-F(clear)^ T • RADSWT. 

CpTT 0(7 

where g is the accelation due to gravity, c p is the heat capacity of air at constant pres- 
sure, RADSWT is the true incident shortwave radiation at the top of the atmosphere (See 
Diagnostic # 48), and 

F(clear)s{y = F(clear)\ w - F(clear)^ sw 


72) OSR Net upward Shortwave flux at the top of the model (Watts/m 2 ) 


rjcp rpNET 

USsK, - r S w,to P 


where top indicates the top of the first model layer used in the shortwave radiation routine. 
In the GEOS-1 GCM, p S w top = 0 mb. 


73) OSRCLR Net upward clearsky Shortwave flux at the top of the model 

(Watts/m 2 ) 

OSRCLR = F(clearsky)sw?t op 

where top indicates the top of the first model layer used in the shortwave radiation routine. 
In the GEOS-1 GCM, psw top = 0 mb. 
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74) CLDMAS Convective Cloud Mass Flux {kg ml sec) 


The cumulative cloud mass flux per unit time, from all convective clouds generated during 
the output interval, is written: 

CLDMAS = rjm B 

where 77 is the entrainment, normalized by the cloud base mass flux, and m B is the cloud 
base mass flux. m B and 77 are defined explicitly in Section 3 . 1 . 1 , the description of the 
convective parameterization. 


75) UAVE Time- Averaged Zonal U-Wind (m/sec) 

The diagnostic UAVE is simply the time- averaged Zonal U-Wind over the NUAVE output 
frequency. This is contrasted to the instantanious Zonal U-Wind which is archived on the 
Prognostic Output data stream. 

UAVE = 


Note, UAVE is computed and stored on the staggered C-grid. 


76) VAVE Time- Averaged Meridional V-Wind (m/ sec) 

The diagnostic VAVE is simply the time-averaged Meridional V-Wind over the NVAVE 
output frequency. This is contrasted to the instantanious Meridional V-Wind which is 
archived on the Prognostic Output data stream. 

VAVE = u(A ,<t>,(T,t) 


Note, VAVE is computed and stored on the staggered C-grid. 


77) TAVE Time- Averaged Temperature {Kelvin) 

The diagnostic TAVE is simply the time-averaged Temperature over the NTAVE output 
frequency. This is contrasted to the instantanious Temperature which is archived on the 
Prognostic Output data stream. 

TAVE = T(A, <f>, <r, t) 


78) QAVE Time- Averaged Specific Humidity {g/kg) 
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The diagnostic QAVE is simply the time- averaged Specific Humidity over the NQAVE out- 
put frequency. This is contrasted to the instantanious Specific Humidity which is archived 
on the Prognostic Output data stream. 

QAVE = <?(A, </>, a, t) 


79) Disabled 


80) PAVE Time-Averaged Surface Pressure - PTOP (mb) 

The diagnostic PAVE is simply the time-averaged Surface Pressure - PTOP over the 
NPAVE output frequency. This is contrasted to the instantanious Surface Pressure - 
PTOP which is archived on the Prognostic Output data stream. 

PAVE = ir(\ r <t>,<j,t) 

= Ps(A,^,cM)-pT 


81) QQAVE Time- Averaged Turbulent Kinetic Energy ( m/sec ) 2 

The diagnostic QQAVE is simply the time-averaged prognostic Turbulent Kinetic Energy 
produced by the GEOS-1 GCM Turbulence parameterization over the NQQAVE out- 
put frequency. This is contrasted to the instantanious Turbulent Kinetic Energy which is 
archived on the Prognostic Output data stream. 

QQAVE = qq(\,(f>,(T,t) 


Note, QQAVE is computed and stored at the “mass-point” locations on the staggered 
C-grid. 


82) SWGCLR Net downward clearsky Shortwave flux at the surface (Watts/m 2 ) 

SWGCLR = F(clearsky)^ nlay+1 

= F ( clearsk y)sW,nh v +l ~ F ( clearsk y)sW,nh v -H 


where nlayTl indicates the lowest model edge-level, or p = p S urf • F (clear sky)SW^ is the 
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downward clearsky Shortwave flux and F(clearsky^gy^ is the upward clearsky Shoitwave 
flux. 


83) AN ALP Analysis Surface Pressure Tendency (mb /sec) 

ANALP is the surface pressure tendency created during the assimilation cycle which is 
used in the Incremental Analysis Updating (IAU) procedure. 

ANALP = (ir A (\,<j),<T,t)- 7r FG '(A>,<M)) / A t Analysis 

where 7 r^ corresponds to the surface pressure after the analysis, tt corresponds to the 
surface pressure first guess, and At Analysis is the time-step between analyses (6 hours). 


84) SDIAG1 User-Defined Surface Diagnostic- 1 

The GEOS GCM provides Users with a built-in mechanism for archiving user-defined di- 
agnostics. The generic diagnostic array QDIAG located in COMMON /DIAG/ , and the 
associated diagnostic counters and pointers located in COMMON /DIAGP/, must be ac- 
ces sable in order to use the user-defined diagnostics (see Section 6). A convenient method 
for incorporating all necessary COMMON files is to include the GEOS GCM vstate.com file 
in the routine which employs the user-defined diagnostics. 

In addition to enabling the user-defined diagnostic (ie., CALL SETDIAG(84)), the User 
must fill the QDIAG array with the desired quantity within the User’s application program 
or within modified GEOS GCM subroutines, as well as increment the diagnostic counter at 
the time when the diagnostic is updated. The QDIAG location index for SDIAG1 and its 
corresponding counter is automatically defined as ISDIAG1 and NSDIAG1, respectively, 
after the diagnostic has been enabled. The syntax for its use is given by 


do j=l,j m 
do i=l,im 

qdiag(i,j,ISDIAGl) = qdiagCi, j ,ISDIAG1) + ... 

enddo 

enddo 

NSDIAG1 = NSDIAG1 + 1 

The diagnostics defined in this manner will automatically be archived by the production 
output programs (see Section 11). 


85) SDIAG2 User-Defined Surface Diagnostic-2 
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The GEOS GCM provides Users with a built-in mechanism for archiving user-defined di- 
agnostics. For a complete description refer to Diagnostic #84. The syntax for using the 
surface SDIAG2 diagnostic is given by 


do j = l , jm 
do i=l,im 

qdiag(i, j ,ISDIAG2) - qdiag(i t j > ISDIAG2) + ... 

enddo 

enddo 

NSDIAG2 = NSDIAG2 + 1 


The diagnostics defined in this manner will automatically be archived by the production 
output programs (see Section 11). 


86) UDIAG1 User-Defined Upper- Air Diagnostic-1 

The GEOS GCM provides Users with a built-in mechanism for archiving user-defined di- 
agnostics. For a complete description refer to Diagnostic #84. The syntax for using the 
upper- air UDIAG1 diagnostic is given by 


do L=l,nlay 
do j=i,jm 
do i=l,im 

qdiag(i, j , IUDIAG1+L-1) = qdiag(i , j ,IUDIAG1+L-1) + ... 

enddo 

enddo 

enddo 

NUDIAG1 = NUDIAG1 + 1 


The diagnostics defined in this manner will automatically be archived by the production 
output programs (see Section 11). 


87) UDIAG2 User-Defined Upper-Air Diagnostic-2 

The GEOS GCM provides Users with a built-in mechanism for archiving user-defined di- 
agnostics. For a complete description refer to Diagnostic #84. The syntax for using the 
upper-air UDIAG2 diagnostic is given by 
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do L=l,nlay 
do j=l , jm 
do i=l,im 

qdiag(i, j , IUDIAG2+L-1) = qdiag(i, j , IUDIAG2+L-1) + ... 

enddo 

enddo 

enddo 

NUDIAG2 = NUDIAG2 + 1 

The diagnostics defined in this manner will automatically be archived by the production 
output programs (see Section 11). 


88) DIABU Total Diabatic Zonal U-Wind Tendency (ml sec /day) 

DIABU is the total time-tendency of the Zonal U-Wind due to Diabatic processes and the 

Analysis forcing. A 

du du du 

DIABU = — + — + — 


Fit 1 m T 'ii.rbii.lp.nce 


F)t 


Note, in Version 1 of the GEOS GCM there are no u-wind tendencies due to cumulus fric- 
tion, (ie. %f tMoist = 0). Also, in forecast simulations ^t Analysis = 0- 


89) DIABV Total Diabatic Meridional V-Wind Tendency (m/ sec/ day) 

DIABV is the total time-tendency of the Meridional V-Wind due to Diabatic processes 
and the Analysis forcing. 


DIABV = 


dv 

dt Moist 



Turbulence 



Analysis 


Note, in Version 1 of the GEOS GCM there are no v-wind tendencies due to cumulus fric- 
tion, (ie. ^ Moist = 0). Also, in forecast simulations §? Analysis = °- 


90) DIABT Total Diabatic Temperature Tendency (deg /day) 

DIABT is the total time-tendency of Temperature due to Diabatic processes and the 
Analysis forcing. 

8T + dT 

Ot MoistProcesses dt ShortwaveRadiation 

dT + dT_ + dT 

dt LongwaveRadiation dt Turbulence dt Analysis 


DIABT = 

+ 
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If we define the time-tendency of Temperature due to Diabatic processes as 


dT dT + dT 

dt Diabatic dt MoistProcesses dt SkortwaveRadiation 

OT +dT 

dt LongwaveRadiation dt Turbulence 

then, since there are no surface pressure changes due to Diabatic processes, we may write 


dT 


p K dirO 


dt Diabatic 7T dt Diabatic 

where 9 = T/p K . Thus, DIABT may be written as 

__ _ . P K ( dirO dirO 

DIABT = — ' 


7T V dt Diabatic dt Analysis 


Note, this diagnostic is only approximate since we have neglected the affect of the analysis 
on the surface pressure, ie. 


<J7T K 


dT _ f_ ( dirO 

dt Analysis 7T \ dt Analysis dt Analysis L P 


+ *% 


1 


This term will be added in subsequent versions of the GEOS GCM. In forecast simulations 
there is no error since . , . =0. 

dt Analysts 


91) DIABQ Total Diabatic Specific Humidity Tendency (g/kg/day) 


DIABQ is the total time-tendency of Specific Humidity due to Diabatic processes and the 
Analysis forcing. 


DIABQ = 


dq dq + dq 

dt MoistProcesses dt Turbulence dt Analysis 


If we define the time-tendency of Specific Humidity due to Diabatic processes as 

dq _ dq ^ dq 

dt Diabatic dt MoistProcesses dt Turbulence 

then, since there are no surface pressure changes due to Diabatic processes, we may write 

dq _ 1 dirq 

dt Diabatic 7T dt Diabatic 

Thus, DIABQ may be written as 

DIABQ = i(^ ) 

7T \ Ot Diabatic Ot Analysis J 
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Note, this diagnostic is only approximate since we have neglected the affect of the analysis 
on the surface pressure, ie. 

dq _ 1 fdirq _ ^ \ 

dt Analysis 7 T \ Ot Analysis dt Analysis / 

This term will be added in subsequent versions of the GEOS GCM. In forecast simulations 
there is no error since f tAnalysis = 0. 

92) Disabled 

93) Disabled 

94) Disabled 

95) Disabled 

96) Disabled 

97) Disabled 

98) Disabled 

99) Disabled 

100) Disabled 

101) Disabled 

102) V1NTUQ Vertically Integrated Moisture Flux (m/sec • g/kg) 

The vertically integrated moisture flux due to the zonal u-wind is obtained by integrating 
uq over the depth of the atmosphere at each model timestep, and dividing by the total mass 
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of the column. 


VINTUQ = 


fsurf U( lP dz 

Isurf P dz 


Using pdz = — ^ 


— -6cr , we have 

9 1 


VINTUQ = / 1 ugdcr 
Jo 


103) VINTVQ Vertically Integrated Moisture Flux (m/sec • g/kg) 


The vertically integrated moisture flux due to the meridional v-wind is obtained by inte- 
grating vq over the depth of the atmosphere at each model timestep, and dividing by the 


total mass of the column. 


VINTVQ = 


f,Zf v< ipdz 

fsurf P dz 


Using pdz = 


— ^ = — -da, we have 

9 9 ’ 


VINTVQ 


[ vqda 
Jo 


104) VINTUT Vertically Integrated Heat Flux (m / sec • deg ) 


The vertically integrated heat flux due to the zonal u-wind is obtained by integrating uT 
over the depth of the atmosphere at each model timestep, and dividing by the total mass 


of the column. 


VINTUT = 


tZf «Tpdz 

Isurf pdz 


Using pdz = —-£- = —jd(T, we have 


VINTUT = f 1 uTdo 
Jo 


105) VINTVT Vertically Integrated Heat Flux (m/sec ‘deg) 
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The vertically integrated heat flux due to the meridional v-wind is obtained by integrating 
vT over the depth of the atmosphere at each model timestep, and dividing by the total 


mass of the column. 


VINTVT = 


Ci 

Cjp dz 


Using p8z — we have 


VINTVT = f 1 vTdcr 
Jo 


106 CLDFRC Total 2-Dimensional Cloud Fracton (0 — 1) 

The 2-dimensional net cloud fraction as seen from the top of the atmosphere is given by 

NLAY 

CLDFRC = 1 - (l - MAX^ ay [CLMOi]) JJ <l-Ci.RO/) 

i=h 

For a complete description of cloud/radiative interactions, see Section 3.2.3. 


107) QINT Total Precipitable Water (gm/cm?) 


The Total Precipitable Water is defined as the vertical integral of the specific humidity, 
given by: 


QINT 



where we have used the hydrostatic relation p8z — 


108) U2M Zonal U-Wind at 2 Meter Depth (m/sec) 


The u-wind at the 2-meter depth is determined from the similarity theory: 


U2M = 


u* u s i ^ 

~k^ m2m W s ~ 


^ 7712 m 

'tftmsi 


Usl 


where , 0 m (2m) is the non-dimensional wind shear at two meters, and the subscript si refers 
to the height of the top of the surface layer. If the roughness height is above two meters, 
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U2M is undefined. 


109) V2M Meridional V-Wind at 2 Meter Depth (m/ sec) 


The v-wind at the 2-meter depth is a determined from the similarity theory: 


V2M = f *^ w . 


V s l Ip. 


m2m 




Vsl 


m st 


where ip m (2m ) is the non-dimensional wind shear at two meters, and the subscript si refers 
to the height of the top of the surface layer. If the roughness height is above two meters, 
V2M is undefined. 


110) T2M Temperature at 2 Meter Depth ( deg K) 

The temperature at the 2-meter depth is a determined from the similarity theory: 

T2M = P K &(i> h2m + t , ) + 0 nrf ) = P«(0 mrf + - 0, urf )) 

* Yh s \ i Yg 

where: 

q* _ { w ' e ') 

u* 

where , 0/i(2m) is the non-dimensional temperature gradient at two meters, 'tpg is the non- 
dimensional temperature gradient in the viscous sublayer, and the subscript si refers to the 
height of the top of the surface layer. If the roughness height is above two meters, T2M is 
undefined. 


Ill) Q2M Specific Humidity at 2 Meter Depth ( g/kg ) 

The specific humidity at the 2-meter depth is determined from the similarity theory: 

Q2M = P*(^(4> hm + %) + q BurJ ) = P«(q, urf + ^ + ^ -(q s! - q, urf )) 

K Vh s i "r tyg 

where: 

(w'q 1 ) 

q* = 

u * 

where iph{ 2m) is the non-dimensional temperature gradient at two meters, ip g is the non- 
dimensional temperature gradient in the viscous sublayer, and the subscript si refers to the 
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height of the top of the surface layer. If the roughness height is above two meters, Q2M is 
undefined. 


112) U10M Zonal U-Wind at 10 Meter Depth {ml sec) 


The u-wind at the 10-meter depth is an interpolation between the surface wind and the 
model lowest level wind using the ratio of the non-dimensional wind shear at the two levels: 


U10M = 


u* u s i 

T pmi0Tn W s 


r m-iom 
V* m s i 




where is the non-dimensional wind shear at ten meters, and the subscript si refers 

to the height of the top of the surface layer. 


113) V10M Meridional V-Wind at 10 Meter Depth (m/sec) 


The v-wind at the 10-meter depth is an interpolation between the surface wind and the 
model lowest level wind using the ratio of the non-dimensional wind shear at the two levels: 




Vsl 


V10M — ^ ^myOm^Y 


$ 


^10 m 




Val 


m sl 


where ^ m (10m) is the non-dimensional wind shear at ten meters, and the subscript si refers 
to the height of the top of the surface layer. 


114) T10M Temperature at 10 Meter Depth {deg K) 

The temperature at the 10-meter depth is an interpolation between the surface potential 
temperature and the model lowest level potential temperature using the ratio of the non- 
dimensional temperature gradient at the two levels: 

T10M = P*(y(ifa 10m + « + <W) = P K (0'urf + - *«“•/)) 

where: 

V * — 

where ^(lOm) is the non-dimensional temperature gradient at two meters, ip g is the non- 
dimensional temperature gradient in the viscous sublayer, and the subscript si refers to the 
height of the top of the surface layer. 
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115) Q10M Specific Humidity at 10 Meter Depth (g/kg) 


The specific humidity at the 10-meter depth is an interpolation between the surface specific 
humidity and the model lowest level specific humidity using the ratio of the non-dimensional 
temperature gradient at the two levels: 

Q10M = + 9W/) = P K (lsur, + - <iw)) 

where: 


where ^(lOm) is the non-dimensional temperature gradient at two meters, if) g is the non- 
dimensional temperature gradient in the viscous sublayer, and the subscript si refers to the 
height of the top of the surface layer. 


116) DTRAIN Convective Cloud Detrainment (kg m/sec) 

The cloud mass flux at the detrainment level per unit time, from all convective clouds 
generated during the output interval, is written: 

DTRAIN = rfa D mB 

where od is the detrainment a level, ra# is the cloud base mass flux, and ij is the entrain- 
ment, defined in Section 3.1.1. 


117) QFILL Filling of negative Specific Humidity (g/kg /day) 

Due to computational errors associated with the numerical scheme used for the advection 
of moisture, negative values of specific humidity may be generated. The specific humidity 
is checked for negative values after every dynamics timestep. If negative values have been 
produced, a filling algorithm is invoked which redistributes moisture from below (see Section 
2.3). Diagnostic QFILL is equal to the net filling needed to eliminate negative specific 
humidity, scaled to a per-day rate: 

QFILL = q%' al 

Q initial 

where 

q n+1 = (wq) n+1 /7r n+1 


118) VINTQANA Precipitable Moisture Adjustment from Analysis ( mm/day ) 
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For the time change in specific humidity due to the analysis increment, § h analysis 1 
vertical integral or total precipitable amount is given by: 


Hop f) q 

P- 


f to P dq dp 7 t f 1 dq 

dz — - I TiT — = - / -^7 da 

is Jsur f at analysis Q g Jo 01 analysis 


VINTQANA = f r 

Jsurf Ot analysis Jsurf 

Defining the specific humidity analysis increment, analysis' 1 as 


dq 


1 (dirq 


Ot analysis 7 T \ dt analysis dt analysis 


dir 


we may write 


VINTQANA =i/‘(^ 


dir , 

-q-frr , . 

analysis Ol analysis , 


119) VINTQFIL Precipitable Moisture Adjustment from Filling (mm /day) 

When advective or other processes create spurious negative specific humidity, total moisture 
is conserved by filling the negative amounts by “borrowing” moisture from below. When 
there is not enough moisture in a column to adequately fill the negative values, the specific 
humidities are set to zero, creating an artificial source of moisture. 

For a change in specific humidity due to filling, A qjuu the vertical integral or total precip- 
itable amount is given by: 

r to P f to P dp tv f 1 

VINTQFIL = / pAqjmdz = - / A qjm— = - / A q S mda 

Jsurf Jsurf 9 9 JO 

A precipitable adjustment rate is defined as the vertically integrated moisture adjustment 
per Dynamics time step, scaled to mm /day. 
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Part III 


GEOS-1 GCM: User’s Guide 

9 Introduction 

In this section we discuss the general philosophy used in the construct of the GEOS GCM 
system, the necessary components of the system, and a basic User’s Guide on how to set-up 
and run experiments. 


10 GEOS Superstructure 


The superstructure of the GEOS GCM is designed so that the entire model may be invoked 
by a simple subroutine call. With each call of the main GCM driver, all model prognostic 
and diagnostic fields are updated by one timestep. A variety of applications may be built 
around the main GCM driver (GCMDRV), which allow for straight forecast simulations, 
assimilation cyles, and dynamic initialization procedures. It should be emphasized that 
what is meant by the GEOS System is simply the library of GCM subroutines which include 
GCMDRV and those that are subsequently called. These subroutines will be denoted as 
below GCMDRV. Those routines which are above GCMDRV include the Main Program, or 
the Application , and any model Output routines. 


10.1 GEOS Applications 

The GEOS Application can be defined as the Main Program from which the GEOS GCM 
and/or the GEOS Utility subroutines are called. All Applications are User-supplied, that 
is, they do not belong to the GEOS GCM compiled library. The purpose of the GEOS 
Application is a pre-determined experiment, or application , which the User wants to perform 
using the GEOS GCM library. 

For example, let us create an Application which performs the following funtions: 


• Read a model restart from Unitl 

• Integrate the model forward in time for NMAX iterations 

• Write the resulting model restart to Unit2 
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This could be accomplished with the following FORTRAN code: 


CALL RESTART (UNIT1,IREAD,...) 

DO N = 1,NMAX 

CALL GCMDRV {Argl, Arg2, ArgS, ...) 

ENDDO 

CALL RESTART (UNIT2,IWRITE,...) 


Here IREAD and IWRITE are read/write 10 flags used by the GEOS Utility Subroutine 
RESTART to read/write GEOS Model restarts, and Argl, Arg2, Arg3 , ... are the the 
arguments to GCMDRV. Both of these argument lists are described fully in the next 
section. 

We can easily modify this application to run forecasts from the same initial state but using 
different GCMDRV model parameters: 


CALL RESTART (UNIT1, IREAD,...) 

DO N = 1,NMAX 

CALL GCMDRV (Arflfl.l, Argl.2, Argl.3 , ...) 
ENDDO 

REWIND UNIT1 

CALL RESTART (UNIT1, IREAD,...) 

DO N = 1,NMAX 

CALL GCMDRV (Arg2.1,Arg2.2,Arg2.3,...) 
ENDDO 

CALL RESTART (UNIT2,IWRITE,...) 


Finally, we could run the model with the same set of GCMDRV model parameters but 
starting from two different initial states: 


CALL RESTART (UNIT!, IREAD,...) 
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DO N = 1,NMAX 

CALL GCMDRV (Arg 1, Arg2 , ...) 

ENDDO 

CALL RESTART (UNIT2,IWRITE,...) 
CALL RESTART (UNIT3,IREAD,...) 

DO N = 1,NMAX 

CALL GCMDRV (Argl, Arg2, Arg3 , ...) 

ENDDO 

CALL RESTART (UNIT4,IWRITE,...) 


These simple examples show the ease and flexibility of creating applications to perform 
experiments with the GEOS GCM. In the next section we describe fully the argument 
lists for subroutines GCMDRV and RESTART, and briefly describe GCMDRVs flow 
structure. 


10.2 Subroutines GCMDRV &; RESTART 

Subroutine GCMDRV is the single interface between User Application routines and the 
actual updating of model prognostic fields. Its flow structure can be separated into three 
essential parts: 


• Read and refresh Boundary Conditions 

• Compute time-tendencies of atmospheric state variables due to phys- 
ical and hydrodynamical processes 

• Integrate the atmospheric state variables forward in time by one 
timestep 


Every call to GCMDRV updates the boundary conditions to the current date and time. In 
Version 1 of the GEOS GCM this is accomplished by a linear interpolation between monthly 
mean data sets. Subsequent versions of the GCM will allow for boundary condition data 
sets with varying frequency intervals. The boundary conditions which are currently updated 
are: 


• Sea Surface Temperature 
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• Ground Wetness 

• Surface Albedo 

• Surface Roughness 

• Sea-Ice Thickness 


In addition, data is provided describing the model topography, topography variance, and 
the Land/Water/Glacier flags. These data are not interpolated to the current date and 
time since they are fixed quantities. Finally, the incident solar radiation at the top of the 
atmosphere is computed for the current date and time. 

At prescribed time intervals, GCMDRV calls separate mini-drivers to update the time- 
tendencies of the prognostic state variables due to particular physical or hydrodynamical 
processes. Each mini-driver has its own timescale of operation. The mini-drivers associated 
with physical and hydrodynamical processes are: 


• Moist Processes (Updated every 10 minutes) 

• Shortwave Radiation (Updated every 3 hours using a normalized in- 
cident radiation, and adjusted every timestep using the true incident 
radiation) 

• Longwave Radiation (Updated every 3 hours) 

• Turbulence (Updated every 30 minutes) 

• Eulerian Hydrodynamics (Updated every timestep) 


The time-tendencies for the processes stated above are always kept in core during the GEOS 
GCM simulation. Thus, all prognostic fields are time-continuous while the time-tendencies 
of the prognostic fields due to the various processes are discontinous at the above stated 
frequencies. These time- tendencies are summed together to provide a total time-tendency 
to be used in the time-integration scheme. Currently, the GEOS GCM can be run with 
either the Matsuno or the Leapfrog time scheme. The Leapfrog scheme has associated with 
it an Asselin (1972) time filter. 
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The calling sequence associated with GCMDRV and a description of the argument list is 
provided below: 

CALL GCMDRV ( KBC, N, NYMD, NHMS, SCHEME, NDT, ALPHA, QQINIT,QBUD ) 


PARAMETER TYPE DESCRIPTION 


KBC 


N 

NYMD 

NHMS 

SCHEME 


NDT 

ALPHA 

QQINIT 


QBUD 


Integer KBC is the unit from which Boundary Conditions are 

read. Every call to GCMDRV will initiate an update of 
the boundary conditions. Boundary conditions are up- 
dated by a linear interpolation between the two months 
kept in memory which strap the current date and time. 
Each time a new mid-month is passed or a new KBC value 
is encountered, two new months are read in. Setting KBC 
equal to zero disables all boundary condition updates. 

Integer N refers to the current time index, and is equal to 1 or 

2. Its value always points to the updated fields and is 
modified by GCMDRV. 

Integer NYMD is the current year/month/day in YYMMDD for- 

mat. Its value is modified by GCMDRV. 

Integer NHMS is the current hour/minute/second in HHMMSS 

format. Its value is modified by GCMDRV. 

Character*4 SCHEME is the variable indicating the time scheme to be 
used for the integration. SCHEME may be set to ‘MATS’ 
for the Matsuno time scheme, or ‘LEAP’ for the Leapfrog 
time scheme. The SCHEME variable may be changed at 
any time during the forecast. 

Integer NDT is the timestep (in seconds) used for the forecast. 

NDT may be positive, zero, or negative, and may be 
changed at any time during the forecast. 

Real Asselin Filter Coefficient, a, to be used with the Leapfrog 

Time Scheme, where q^ tlt = q n ( 1 - a) + a ^ 3a±i±gii=i ^ . 

ALPHA may be changed at any time during the forecast. 

Logical QQINIT is a Logical flag parameter used by the Turbu- 

lence parameterization scheme. QQINIT should be set to 
TRUE if the GEOS GCM is starting from an initial con- 
dition which does not contain Turbulent Kinetic Energy 
as a prognostic state variable. After one successful call to 
GCMDRV, QQINIT should be set to FALSE. 

Logical Logical flag to compute budget diagnostics. 
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The GEOS GCM Utility subroutine RESTART may be used within User-Defined Applica- 
tions to read and/or write a GEOS GCM Restart dataset. The calling sequence associated 
with RESTART and a description of the argument list is provided below: 


CALL RESTART ( KU,IOFLAG,NYMD,NHMS,N,NYMDO,NHMSO,JOB,NTYPE ) 


PARAMETER 

TYPE 

KU 

Integer 

IOFLAG 

Integer 

NYMD 

Integer 

NHMS 

Integer 

N 

Integer 

NYMDO 

Integer 

NHMSO 

Integer 

JOB 

Character*8 

NTYPE 

Integer 


DESCRIPTION 


Unit number from which to read or to write the GEOS 


10 Flag which is equal to: 1 to Read a GEOS Restart, 2 

to Write a GEOS Restart. . , 

NYMD is the current year/month/day m YYMMDD for- 


NHMS is the current hour/minute/second in HHMMSS 

format. . . , . . , , 1 

N refers to the current time index, and is equal to 1 or 

2. Its value always points to the most recently updated 

NYMDO is the beginning year/month/ day of the current 
experiment in YYMMDD format. NYMDO is used in cur- 
rent Production Applications to control Output frequen-’ 

Gigs . 

NHMSO is the beginning hour /minute/second of the cur- 
rent experiment in HHMMSS format. NYMDO is used 
in current Production Applications to control Output fre- 
quencies. 

An Experiment Identifier (usually a name or experiment 
number). 

NTYPE is the GEOS Restart Identification index. 
NTYPE = -1 refers to a “First Guess” restart which was 
used for an analysis. NTYPE = 0 refers to the “After 
Analysis” or straight forecast restart. NTYPE = 1 refers 
to a restart from an IAU assimilation. 


It should be noted that if IOFLAG = 1, the RESTART parameters NYMD, NHMS, N, 
NYMDO, NHMSO, JOB, and NTYPE are read from the GEOS Restart dataset and passed 
back to the calling program. If IOFLAG - 2, the parameters NYMD, NHMS, N, NYMDO, 
NHMSO, JOB, and NTYPE passed into subroutine RESTART will be written out to the 
GEOS Restart. 
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11 GEOS Production 


Although the GEOS GCM system is designed to work with a User supplied “Applications” 
program, several application programs have been developed by the DAO which are used 
in standard production forecast simulations and analyses, and are available as templates. 
These applications provide examples of convenient ways to control the input and output 
streams, the frequency and selection of diagnostic and prognostic fields for output, and 
the simulation start and end times. The application routine examples reside on the SGI 
workstation hera.gsfc.nasa.gov , in the directory /production/geos-das/geos 1.1 /applications. 
Some of them are described here in detail: 

gcmprod.f - A simple application routine designed to read and follow namelist-supplied 
directions to simulate a specific period and direct all “activated” diagnostics 
to an output stream. The application reads one namelist to establish the 
initial and final time and date of the simulation, the frequency of the out- 
put stream, and other run control parameters such as timestep and asselin 
filter strength. Another namelist is read in which logical flags determine the 
selection of diagnostics to be activated and directed to the output stream. 

gcmprodn.f - An application routine which is designed to handle a somewhat more complex 
output stream for diagnostic and prognostic fields. The namelist controlled 
interface for the diagnostics involves the specification of an HHMMSS for- 
mat variable for each activated diagnostic separately. Diagnostics may be 
averaged for any time period up to the length of the simulation, and the av- 
eraging period may be centered about the current time or may extend from 
the previous time to the present. The application determines how many dif- 
ferent averaging periods and types are present, and directs the output to the 
appropriate number of separate output streams. Separate output routines 
are called to properly direct the diagnostic output. Control over the begin 
and end times of the simulation are also governed by namelist parameters. 

amip.f - An experiment-specific application routine that was designed to handle the 

control of an extended multi-year simulation with complex diagnostic output 
requirements. The application was written to direct diagnostic output to 
several output data streams each with a pre-specified set of diagnostics. 
Separate output routines are called for each data stream, providing them 
with the required list of diagnostic quantities. 


In addition to flexible application routines, the GEOS GCM system also incorporates user 
defined output routines that can be called by the application. The output routines developed 
at the DAO write a combination of prognostic and/or diagnostic quantities, on model sigma 
levels or interpolated to pressure levels, in the Phoenix Output Format (see Section 11.2). 
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Just like the application routines, several output routines have been developed by the DAO 
and are available as examples. Some of them are described here: 


diagprs.f 


diagsig.f 


glaoutc.f 


sigoutc.f 


progprs.f 


progsig.f 


An output routine which accepts as arguments a specific list of diagnostic 
quantities and interpolates them to an argument specified list of pressure lev- 
els. The routine uses the model provided GETDIAG subprogram to obtain 
the proper time-averaged quantities. The code is written to distinguish be- 
tween different quantities and attempts to use an appropriate interpolation 
algorithm. Output is written in Phoenix Format. 

Designed to accept an argument- specified list of diagnostics, and direct out- 
put in Phoenix Format to one output stream. The GETDIAG routine is 
used to retrieve diagnostic quantities in the proper units from the model. 

When this output routine is called, all prognostic quantities and all activated 
diagnostic quantities are interpolated to an argument-list specified set of 
pressure levels. All output is directed to a Phoenix Format data stream. 

All model prognostic and activated diagnostic quantities are directed to a 
single sigma- level Phoenix Output Format data set. 

Designed to direct model prognostic quantities to a single output stream. 
Quantities are interpolated to a list of pressure levels. 

Model prognostic quantities are simply directed to a single sigma-level 
Phoenix format data set. 
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11.1 Production Namelists 


Two general namelists have been designed to setup and control the Production runs of the 
GEOS GCM. The namelists include information about the beginning date and time of the 
forecast, the ending date and time of the forecast, the duration of each job segment, and 
output characteristics. 

The first namelist is called INPUT, and contains the parameters to control the beginning 
and ending times of each experiment, as well as some output information. The following 
INPUT namelist and subsequent description can be used as is or as a template for Users to 
build upon within their own Applications: 


&INPUT 








JOB 

= 

’ 1234’ , 






XLAB 

= 

’E1234 GEOS 

-GCM 

CLIMATE 

SIMULATION EXPE’ 

) 



’RIMENT 





9 

NYMDB 

= 

850101 , 






NHMSB 

= 

000000 , 






NYMDE 

= 

860101 , 






NHHSE 

= 

000000 , 






NDT 

= 

225 






LEAP 

= 

.TRUE. , 






MATS 

= 

.FALSE. , 






QQINIT 

- 

.FALSE. , 






ALPHA 

= 

0.05 






NDSEG 

= 

2400000 , 






NDOUT 

= 

060000 , 






PLEVS 

= 

1000.0, 

950.0 

, 900.0, 

850.0, 

800.0, 




700.0, 

600.0 

, 500.0, 

400.0, 

300.0, 

250.0, 



200.0, 

150.0 

, 100.0, 

70.0, 

50.0, 

30.0, 


&END 
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PARAMETER TYPE 


DESCRIPTION 


JOB Character*8 Experiment Identifier (Usually a name or experiment 

number) 

XLAB Character *80 User Assigned Description of Current Experiment 

NYMDB Integer Beginning Date of Experiment in YYMMDD format 

NHMSB Integer Beginning Time of Experiment in HHMMSS format 

NYMDE Integer Ending Date of Experiment in YYMMDD format 

NHMSE Integer Ending Time of Experiment in HHMMSS format 

NDT Integer Time step (in seconds) of the hydrodynamics 

LEAP Logical Logical Flag to use the Leapfrog Time Scheme 

MATS Logical Logical Flag to use the Matsuno Time Scheme 

QQINIT Logical QQINIT is a Logical flag parameter used by the Turbu- 

lence parameterization scheme. QQINIT should only be 
set to TRUE if the GEOS GCM is starting from an initial 
condition which does not contain Turbulent Kinetic En- 
ergy as a prognostic state variable. After one successful 
call to GCMDRV, QQINIT should be set to FALSE. 
ALPHA Real Asselin Filter Coefficient, a, to be used with the Leapfrog 

Time Scheme, where q^ llt = q n ( 1 - «) + a (“ n+1 

NDSEG Integer Current job segment time interval in HHMMSS for- 

mat. The total time of a GEOS simulation is given as 
(NYMDE, NHMSE) - (NYMDB, NHMSB). This total time 
is divided into a number of smaller job segments, which 
run in succession on the computer. The job segment size 
is determined by the User, depending on the desired size 
of the Output stream and/or the CPU time restrictions 
within the computer QUEUE. Note, HH may be any num- 
ber of hours, and not limited to a 24 hour day (ie., 5 day 
job segments would be defined as NDSEG = 1200000). 
NDOUT Integer Output frequency for Prognostic fields in HHMMSS for- 

PLEVS Real Levels used for Sigma to Pressure Interpolation to create 

Pressure Level Output. (Maximum number of levels - 
50) 

When the GEOS GCM is run during for a particular experiment, no diagnostic quantities 
will be produced or saved unless the User has explicitly enabled the desired diagnostics fiom 
the Applications program. Using the DAO supplied applications directly or as templates, 
a convenient NAMELIST method has been adopted for diagnostic control. Diagnostic 
quantities which are of interest to the User may be enabled by the application by including 
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the Diagnostic Name from the Diagnostic Menu (Section 7) in the namelist INDIAG. For 
those Users employing the gcmprod.f application described in Section 11, simply preface the 
Diagnostic Name with the letter L. This variable will then be treated as a logical parameter. 
Those parameters set to TRUE will be enabled for the GEOS GCM simulation. An example 
of the INDIAG namelist using the Logical variable syhtax is shown below: 


&INDIAG 



LUFLUX 

= 

.TRUE. , 

LVFLUX 

= 

.TRUE. , 

LHFLUX 

= 

.TRUE. , 

LMOISTT 

= 

.TRUE. , 

LMOISTQ 

= 

.TRUE. , 

LRADLW 

= 

.TRUE. , 

LRADSW 

= 

.TRUE. , 

LPREACC 

= 

.TRUE. , 

LPRECON 

= 

. TRUE . , 

LQLR 

= 

.FALSE. 

LEVAP 

= 

.FALSE. 

LU2M 

= 

.FALSE. 

LV2M 

= 

.FALSE. 

LT2M 

= 

.FALSE. 

LQ2M 

= 

.FALSE. 


&END 


In this example, the diagnostic quantities UFLUX, VFLUX, HFLUX, MOISTT, MOISTQ, 
RADLW, RADSW, PREACC, and PRECON will all be enabled during the simula- 
tion. The remaining diagnostics OLR, EVAP, U2M, V2M, T2M, and Q2m will be 
disabled during the run. It should be noted that explicitly defining the logical diagnostic 
parameters to FALSE in the INDIAG namelist is not necessary since this is the default. 
Using this method in conjunction with the application gcmprod.f will produce diagnostic 
quantities whose output frequency is governed by the INPUT namelist parameter NDOUT. 

For those Users interested in a more sophisticated diagnostic output frequency algorithm, 
each diagnostic may be written out at its own unique frequency by using the DAO applica- 
tion gcmprodn.f In this case, a particular diagnostic is turned on, or enabled, by specifying 
a non-zero diagnostic frequency, in HHMMSS format, for that diagnostic. The diagnostic 
frequency name is the same as the diagnostic name (from the Diagnostic Menu, Section 
7) preceded by the letter N (eg, NUFLUX = 030000). In this example, the diagnostic 
UFLUX will be accumulated and written out every 3 hours. Every diagnostics has its own 
frequency, thus certain diagnostics may be written out more frequently than others. A 
separate file is defined for each unique diagnostic frequency. Thus, if the user turns on a 
total of 24 diagnostics, 12 with a 6 hour frequency, 10 with a 3 hour frequency, and 2 with a 
half-hour frequency, a total of 3 diagnostic datasets will be created. Finally, diagnostics are 
accumulated and written out with the time-stamp at the end of the accumulation period. 
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A time-stamp which occurs at the middle of the accumlating period may be obtained by 
making the diagnostic frequency negative. An example of the INDIAG namelist using the 


frequency 

variable sy 

&INDIAG 

NUFLUX 

= 030000, 

NVFLUX 

= 030000, 

NHFLUX 

= 030000, 

NMOISTT 

=-060000, 

NMOISTQ 

=-060000, 

NRADLW 

=-060000, 

NRADSW 

=-060000, 

NPREACC 

= 060000, 

NPRECON 

= 060000, 

NOLR 

= 060000, 

NEVAP 

= 060000, 

NU2M 

= 010000, 

NV2M 

= 010000, 

NT2M 

= 010000, 

NQ2M 

&END 

= 010000, 


In this example, four output data streams will be produced. Diagnostics U2M, V2M, 
T2M, and Q2M will be averaged over a one-hour frequency and written out with a time- 
stamp at the end of each hour. Diagnostics UFLUX, VFLUX, and HFLUX will be 
averaged over a three-hour time-period. Diagnostics PREACC, PRECON, OLR, and 
EVAP will be averaged over a six-hour period and written out with a time-stamp at the 
end of the averaging period. Diagnostics MOISTT, MOISTQ, RADLW, and RADSW 
will also be averaged over a six-hour period but will be written out with a time-stamp 
associated with the middle of the averaging period. All other diagnostics are disabled. 
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11.2 Phoenix Format 


The Phoenix Format used in the archiving of the GEOS AMIP and the GEOS-DAS 5- Year 
Re- analysis employs a self- documenting output stream which describes each field and its 
attributes for each time period on the output dataset. At each output time period, Header 
information is given which describes the horizontal and vertical dimensions of the output, 
the time period of the output, and variable descriptions of the data type. The number and 
type of fields are User Defined. 

The Phoenix format allows for writing data in three different ways (or configurations). The 
first, configuration #1, is simply a sequence of distinct 2-dimensional fields (such as sea- 
level pressure, topography, or selected individual levels of specific quantities such as wind or 
moisture). The second method, configuration #2, is 3-dimensional data which are ordered 
by quantity for each level. It is assumed that all quantities written in this manner are 
defined for all levels. The third, configuration #3, is 3-dimensional data which are ordered 
by level for each quantity. Here, the number of levels associated which each quantity may 
vary. The specific configurations used within any dataset are determined by the User. The 
following sample code may be used to read any Phoenix Format Output stream: 


C *************************************************************************** 
C **** GODDARD LABORATORY FOR ATMOSPHERES **** 
C *************************************************************************** 


PARAMETER ( IDIM = 144 ) ! IDIM=144 for 2.5 degree Longitude 

PARAMETER ( JDIM = 091 ) * JDIM=091 for 2.0 degree Latitude 

DIMENSION Q( IDIM, JDIM ) 

CHARACTER+8 XLABEL(IO) , JOB, NAME 
DIMENSION SIGE (50) , NDLEV(IOO), ZLEV(50) 

DIMENSION IDUM(IOO) , RDUM(IOO) 

CHARACTER+8 NAMES(lOO), NAMEU(IOO), NAMED ( 100) , CDUM(100) 

CHARACTER+40 DESCS(IOO), DESCU(IOO), DESCD(IOO), DDUM(100) 


C *************************************************************************** 
C **** READ HEADER **** 
C *************************************************************************** 


DATA KU /10/ 

READ (KU , END =5 00) JOB, NYMD , NHMS, NYMD0 , NHMS0, 

XLABEL, IM, JM, 

NSFLD , NUFLD , NDFLD , 

PT0P, NULEV, (ZLEV(K) ,K S 1 ,NULEV) , 

NLAY, (SIGE(K) ,K=1 ,NLAY+1) , 

NDUM , (DDUM(N) , IDUM(N) , RDUM ( N ) , CDUM ( N ) , N= 1 , NDUM ) 


READ(KU) ( NAMES (N) , DESCS(N), 


N=l, NSFLD ), 
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( NAMEU(N) , DESCU(N) , N=1 ,NUFLD ), 

( NAMED (N) , DESCD(N) , NDLEV(N) , N=1 ,NDFLD ) 


Q *************************************************************************** 

c **** READ 2-DIMENSIONAL FIELDS **** 

C ********* ********************************************************* ********* 


DO N=1 ,NSFLD 

READ(KU) ZLEVID,NAME,Q 

ENDDO 


C *************************************************************************** 
c **** READ 3 -DIMENSIONAL FIELDS (ORDERED ALL FIELDS FOR EACH LEVEL) **** 

C *************************************************************************** 

DO L= 1 , NULEV 

DO N=1 ,NUFLD 

READ(KU) ZLEVID,NAME,Q 

ENDDO 

ENDDO 

C ************************ ********* ******************* *********************** 

C **** READ 3-DIMENSIONAL FIELDS (ORDERED ALL LEVELS FOR EACH FIELD) **** 

C ***************************************************** ************ ********** 


DO N=1 ,NDFLD 

DO L=1 ,NDLEV(N) 

READ(KU) ZLEVID ,NAME,Q 

ENDDO 

ENDDO 

500 CONTINUE 
STOP 
END 
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PARAMETER 

TYPE 

DESCRIPTION 

JOB 

Character*8 

Experiment Identifier (Usually a name or experiment number) 

XLAB 

Character*80 

User Assigned Description of Current Experiment 

NYMD 

Integer 

Current Date of Experiment in YYMMDD format 

NHMS 

Integer 

Current Time of Experiment in HHMMSS format 

NYMDO 

Integer 

Beginning Date of Experiment in YYMMDD format 

NHMSO 

Integer 

Beginning Time of Experiment in HHMMSS format 

IM 

Integer 

Longitudinal dimension of data 

JM 

Integer 

Latitudinal dimension of data 

NLAY 

Integer 

Number of GCM vertical levels used in Experiment 

SIGE 

Real 

GCM Sigma-Edge values 

PTOP 

Real 

Model Top Pressure used in Experiment 

ZLEV 

Real 

Upper- Air Level Values archived on data (Pressure or Sigma 
Level) 

NSFLD 

Integer 

Number of 2-Dimensional Fields written in Config. #1 

NUFLD 

Integer 

Number of 3-Dimensional Fields ordered by Level written in 
Config. #2 

NDFLD 

Integer 

Number of 3-Dimensional Fields ordered by Quantity written 
in Config. #3 

NAMES 

Character*8 

Names of 2-Dimensional Fields defined by NSFLD 

DESCS 

Character*40 

Description of 2-Dimensional Fields defined by NSFLD 

NAMEU 

Character*8 

Names of 3-Dimensional Fields ordered by Level defined by 
NUFLD 

DESCU 

Character*40 

Description of 3-Dimensional Fields ordered by Level defined 
by NUFLD 

NULEV 

Integer 

Number of Upper- Air Levels associated with all fields in Config. 
#2 

Names of 3-Dimensional Fields ordered by Quantity defined by 
NDFLD 

NAMED 

Character*8 

DESCD 

Character*40 

Description of 3-Dimensional Fields ordered by Quantity de- 
fined by NDFLD 

NDLEV 

Integer 

Number of Upper- Air Levels associated with each field in Con- 

fig- #3 

NDUM 

Integer 

Number of User-Defined Header variables. NDUM, DDUM, 
IDUM, RDUM, and CDUM may be used by Users to add 
REAL, INTEGER, or CHARACTER information in the Header 
concerning the User’s experiment. 

DDUM 

Character*40 

User-Defined Description of Header variables 

IDUM 

Integer 

User-Defined Integer Header variables 

RDUM 

Real 

User-Defined Real Header variables 

CDUM 

Character*8 

User- Defined Character Header variables 

PLEVS 

Real 

Levels used for Sigma to Pressure Interpolation to create Pres- 
sure Level Output. (Maximum number of levels = 50) 

ZLEVID 

Real 

Level Indentifier (Pressure or Sigma Level) of Quantity being 
read 

NAME 

Character*8 

Character Name of Quantity being read 

Q 

Real 

Two-Dimensional Quantity being read 
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It should be noted that the Datasets produced for the GEOS-1 DAS 5-year re-analysis and 
the GEOS-1 GCM AMIP simulation have the following characteristics and User-Defined 
Header information: 


• The horizontal resolution for all data is constant for all time periods. 

• All horizontal data, q(i, j),is written with i = 1 corresponding to Longitude - 
-180 (International Dateline), and j = 1 corresponding to Latitude = -90 
(South Pole). 

t The data written in Configuration # 1 are generally surface data (such as 
Sea-Level Pressure, Surface Topography, Land/ Water Masks, etc.). 

• The data written in Configuration # 2 are generally Upper- Air prognostic 
quantities (such as wind, temperature, moisture, etc.). 

• The data written in Configuration # 3 are generally GCM Model Diagnostic 
quantities (such as Diabatic Heating). Some Diagnostics are only defined at 
one level (such as Precipitation, Outgoing Longwave Radiation, etc.). 


The User-Defined Header information has been defined as: 


NDUM 

= 2 

DDUM(l) 

= ’OUTPUT RECORD TYPE (INTEGER)’ 

IDUM(l) 

= -1 FIRST GUESS 

0 AFTER ANALYSIS (or STRAIGHT FORECAST) 

1 INITIALIZED ANALYSIS 

RDUM(l) 

= 0.100000E+16 

CDUM(l) 

= ’LOG8R’ 

DDUM(2) 

= ’UNDEFINED FILL VALUE (REAL)’ 

IDUM(2) 

= 0 

RDUM(2) 

= 0.100000E+16 

CDUM(2) 

= ’UNDEF’ 
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11.3 Geopotential Heights on Sigma Surfaces 

In the GEOS-1 GCM, the pressure gradient computed on sigma surfaces is derived by 
removing the large-scale component due to a mean potential temperature and surface pres- 
sure. Therefore, the perturbation geopotential heights have been archived on the GEOS-1 
GCM sigma level datasets. 

The pressure gradient on a constant pressure surface is related to the gradient on a constant 
sigma surface by 


V p $ = V a <I> + c p 6V a p K , 
while the hydrostatic equation is given as 


< 94 > 

dp K 


— c p 0 . 


Defining $ = $ + and 0 = 0 + 0', we may write 


d ($ + $') 

dp K 



(25) 


(26) 


(27) 


or 


dp K 


— c p 0 


$ = -c p 0{p K -Po) 


dp K 


- c p 0 ' . 


Thus, we may express the pressure gradient on sigma surfaces as 


(28) 

(29) 


+ c p 0V a p K 


v ($ + $') +c p (e + »’) vp K 
v<E + v$' + c„evp K + c p e'v P K 

-c p 0Vp K + V$' + c p 0Vp K + cpB'Vp* 
V$' + c p 0'Vp K 


( 30 ) 
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where 


= $t- c p 0 (pg - pf) 

0'i = 6i-0 ; 0 = 280/1000 K 

po = 1000 . (31) 
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11.4 Interpolation from Sigma to Pressure Surfaces 


While the GEOS GCM uses the sigma coordinate in the vertical for the discritization of 
the primitive equations, the current analysis system requires “first guess” fields to be on 
prescribed pressure surfaces. In addition, it is often desirable to process model output on 
pressure surfaces for scientific studies and comparisons to observations. The method used 
for the sigma coordinate to pressure coordinate interpolation is consistant with the form of 
the hydrostatic equation used in the GEOS GCM. 

The vertical grid and index notation for sigma levels used in the GEOS GCM is given by 


f. A -R. t- 

q* 

ie 


$ P* 9 q_ 


Figure 6: Vertical placement and index notation for sigma levels in the GEOS GCM 
The form of the hydrostatic equation used in the GEOS GCM is given by 


9* __ „ 

dP n - C P 8 


(32) 


where $ is the geopotential height, and 6 = T/P K . Similarly, for an arbitrary quantity, q , 
we may define 


dq 


dpn 


— c 


(33) 


where c is a constant. Thus, we may write 


dq_ 

dP K 


= c 


qi-\ — qi qe-i — q * 


DK D/C p/t 

r t - 1 l t~\ 


pK 


(34) 


or 
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q* = qi - 1 


(35) 


1- 


DK DK 

*1-1 ~ J * 

DK 

Pf-1 r £ . 


DK DK 

-*1-1 

J?-i - P( 


This method is used for the winds, temperature, and relative humidity. Specific humidity 
at pressure levels is obtained from the previously interpolated temperature and relative 
humidity. For target pressures above the first model level, a linear extrapolation in P K 
from levels 1 and 2 is used. For target pressures between the surface pressure and the 
lowest model level, the quantity is assumed constant and is equal to its value in the lowest 
model level. For target pressures greater than the surface pressure (below topography), the 
quantity is undefined. It should also be noted that all GEOS GCM diagnostic quantities 
on pressure surfaces have been similarly linearly interpolated in P K . 

The method used in the interpolation of geopotential heights, in particular, must be inti- 
mately linked to the finite- differencing scheme used in the background model in order to 
minimize computational error and prevent spurious height gradients (or geostrophic imbal- 
ances). The technique described here, formulated at UCLA by Arakawa and Suarez (per- 
sonal communication), uses two separate interpolation algorithms depending on whether 
you are above or below the nearest sigma-level interface, fe, and has the following proper- 
ties: 


$ P K 6 


© 

P K 0 

ie. 


<f> P K 0 

1 


Figure 7: Index notation for vertical interpolation of Geopotential Heights 


1) 

2 ) 

3 ) 

4 ) 

5 ) 


$Pl = §1 

$p 2 = 
y p — 

VpWi = V,*, + CpOtVrP? 

V p $p 2 = + Cp0z-.iV, 


for Pi = Pi 
for P 2 = Pi - 1 
for Pi = P 2 = Pie 
for Pi = Pi 
Pi - 1 for P 2 = Pi - 1 
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Thus, the interpolation not only preserves the heights on sigma surfaces whose target pres- 
sure coincides with the sigma level pressure, but also the horizontal gradient of the height 
on pressure surfaces is defined accurately as the horizontal gradient of the sigma level height 
plus the term corresponding to the horizontal gradient of the surface pressure. 

We begin by writing the hydrostatic equation in the form 


dP K 


— Cp9 . 


Following Arakawa and Suarez (1983), we may write 


(36) 


Ote 

Thus, for the two regions in 
1) For Pe > Pi > Pie 

= <t>( + c p 0 te (P? - 

(PPe ~ P?-l)9t- 

-PLih 

■i + {Pi - p&M 

(37) 

(38) 

(Pf- 

Figure 7, we have: 

-PiU) 


$Pi = $i + c p f 
JP 

p t 

&ldp K 

1 

(39) 

2) For Pi e > P 2 > Pi- 1 

$P 2 = $i-i + c p 

p 2 

e 2 dP K 

(40) 


where 


®i(P K ) 

02(P*) 


(Pf -P")9 fe + (P"- 

(p?'-p*)e t . 1 + (p«-p e *_ 1 ] 

(P£ - PL i) 


(41) 

(42) 


Here 0i(P K ) and 0 2 (jP' c ) are the interpolated 0 values from levels i and t — 1 to the pressure 
P. Using (42) in (39), we may define 
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( 43 ) 



r F > (Pr-P^he + jP'-H) 6 ' ^* 

JPi (H - He) 


h = 


he {pm - Pi] - \\Pf ~ *T 2 ]) 

+e t {\[ph - pf] - pm - pn) 


Pf-P\ 

Pt-Pu 


he {p* - \[P? + Pli) 
+0e Qt Pi + Pi] - He) 


(44) 


(45) 


1 

2 


P/-P, K 
P* P(e 


he (Pi -P?)~ »t (m - Pi - Pi) 


(46) 


Similary, we may define 


or 



[ P > (PZ-P^Ol-l+jP*- P?-l)he (!r , 

Jpe-i ( P £ e ~ P i-l) 


(47) 


h = 


9 t - 1 {pm - pl ii - \(p f - p t- i 2 i) 

+he { \[P f - Pl- 1 2 ] - PI m - Pl- ll) 


(48) 
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Pj-Pj-I 

Pf-Pf , 


»t-l [Pie ~ 7^2 + PL l] 


( 5 [ PS + Jfc,] - PL 1 


1 ( P?- p t-y 

2 W 4 -P,li 


-1 (2P£ - Pi - PLi) + 4 (P? - P£-i)' 


Therefore, for the interpolation of heights within regions (1) and (2) from Figure 7, we may 
write: 

1) For Pi > P y > Pi e 


^ ^ c P / P? ~ Pi 

<Fp, = — A 

1 2 l P" - P« 


(Pi - Pi) - <k (iPte - Pe - Pf)] 


2) For P ee > P 2 > Pe - 1 


<Fp, = $f_i + 


Pi -Pi 


2 \p;; - Pi 


‘rj K-i (2 Pie - Pi - PLi) + 4 (Pi - PLi)} 
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11.5 Sea-Level Pressure 


In addition to the archiving of surface pressure, which is a prognostic variable in the GEOS- 
1 GCM, on all GEOS Production output streams, the sea-level pressure is also computed 
and archived. Also, during the GEOS data assimilation cycle, sea-level pressure rather than 
surface pressure is analyzed. The method used to compute sea-level pressure from surface 
pressure involves integrating the hydrostatic equation from sea-level to the surface, and an 
extrapolation of the temperature profile below topography. 

The hydrostatic equation may be written as: 


dp pg 

dz~ P9 ~ RT ' 

Integrating (53) from sea-level, s/, to the surface, s, we may write 


(53) 



r s d$ 

h ST 


(54) 


where we have used = gdz. Assuming a mean temperature, T, between the sea level 
and surface, (54) may be written as 


or 


In(^) 

Psl 


RT 


Psl 


= p s e 


$ S /RT 


(55) 

(56) 


The mean temperature T is defined as the average of the surface and sea- level tempeiatuies. 


T = \(T a + T sl ) 


(57) 


where T s i is obtained by assuming the temperature follows the moist adiabatic lapse rate 


dT_ 

dz 


= -/* 


T s i = T S + 




(58) 


where (3 — 6.5° The surface temperature used for the sea- level pressure calculation is 
obtained by averaging the potential temperature, 0, in the lowest 100 mb of the model. 
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This allows variable resolution changes in the PBL while maintaining consistent sea-level 
pressure values between resolutions. We have 


T a - 0p« ; 0 = T/p K 


where 


such that 


£ = El 


Y A a i = 

/L-J 

t 


Ei &Pe 

7T 


100 mb 

IT 


( 59 ) 

( 60 ) 

( 61 ) 


Here Jr = p s - p top - 
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12 GEOS Unix Script 


The GEOS GCM is currently being run on NASA/GSFC’s Cray C-90 computer system 
under the Cray UNICOS environment. There are five essential components to the basic 
script which are needed in order to perform experiments with the GEOS GCM: 


• Create GEOS GCM NAMELIST for model parameters and diagnos- 
tics 

• Create DIAGSIZE Subroutine for enabled diagnostics 

• Remote copy frozen GEOS Libray, Application and Output Routines 
from geos-das@dao 

• Compile User Application and Output Routines 

• Load and Run GEOS simulation 


Each of these elements are briefly discussed in the following subsections. 


12.1 Creating the GEOS GCM Namelist 

The following example shows the creation of a typical GCM namelist to be used during a 
GEOS simulation: 


################################################################## 

###### 

###### Setup GCM Namelist 

###### 

################################################################## 


/bin/rm gcmrun .namelist 
cat > gcmrun. namelist « > E0F’ 


&INPUT 



JOB 

= 

> 1234’, 

XLAB 

= 

, E1234 GEOS-GCM CLIMATE SIMULATION EXPE 



’RIMENT 

NYMDB 

= 

850101 , 

NHMSB 

= 

000000 , 

NYMDE 

= 

860101 , 

NHMSE 

= 

000000 , 

NDT 

= 

225 

LEAP 

= 

.TRUE. , 

MATS 

= 

.FALSE. , 


90 


(TZ. . 



ALPHA =0.05 
NDSEG = 2400000 , 

NDOUT = 060000 , 


PLEVS 

= 

1000.0, 

950.0 

, 900.0, 

850.0, 

800.0 



700.0, 

600.0 

, 500.0, 

400.0, 

300.0 



200.0, 

150.0 

, 100.0, 

70.0, 

50.0 


&END 


&INDIAG 

NUFLUX 

= 030000, 

NVFLUX 

= 030000, 

NHFLUX 

= 030000, 

NM0ISTT 

=-060000, 

NM0ISTQ 

=-060000, 

NRADLW 

=-060000, 

NRADSW 

=-060000, 

NPREACC 

= 060000, 

NPREC0N 

= 060000, 

N0LR 

= 060000, 

NEVAP 

= 060000, 

NU2M 

= 010000, 

NV2M 

= 010000, 

NT2M 

= 010000, 

NQ2M 

= 010000, 

&END 

EOF' 



250.0, 

30.0, 20.0 


For a complete description of GCM Namelist parameters, see Section 11.1. 

12.2 Creating th% DIAGSIZE Subroutine 

Using the DIAGSIZE subroutine is a convenient method to allocate just enough perma- 
nent memory during a run to be used for diagnostic storage. All GEOS GCM diagnostic 
quantities are stored in the single diagnostic array QDIAG which is located in the DIAG 
COMMON. In order to minimize the maximum memory requirement used by the model, 
the GEOS GCM libraries are compiled with room for only one horizontal diagnostic array. 
All production DAO applications contain a call to subroutine DIAGSIZE, which expands 
the size of the DIAG COMMON to include the total number of surface and upper-air di- 
agnostics the User has enabled. For the diagnostics enabled in the preceeding section, the 
appropriate DIAGSIZE subroutine should be defined as: 

################################################################## 

###### 

###### Create DIAGSIZE Routine to control Total Diagnostics 
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###### 

########################################################### ####### 


/bin/rm diagsize.f 
cat > diagsize.f CC’EOF’ 
subroutine diagsize 
parameter ( imr = 144 ) 

parameter ( jmr = 90 ) 

parameter ( nlayr = 20 ) 

parameter ( ndiags = 11 ) 

parameter ( ndiagu = 4 ) 

common /diag/ ndiagmx,qdiag(imr , jmr+1 ,ndiags+ndiagu*nlayr) 

ndiagmx = ndiags + ndiagu*nlayr 

return 

end 

’EOF* 


We see that the User needs to simply assign the appropriate values for the number of 
surface diagnostics (NDIAGS) and the number of upper-air diagnostics (NDIAGU) which 
are enabled in the User’s namelist. 


12.3 Remote copy frozen GEOS Library 

In order to gain access to the GEOS GCM, the User must remote copy the frozen production 
system from the DAO production workstation, gcos-d(is@dcio. In addition to the GEOS 
GCM library, the User may also remote copy production application programs and output 
routines. The COMMON files associated with these programs must also be copied to the 
User’s local directory. Finally, the GEOS GCM boundary conditions and othei GEOS 
input datasets must be obtained from the production workstation. An example of the 
syntax necessary for these remote copies is as follows: 


################################################################## 

###### 

###### Get Frozen GCM, Applications, and Output Routines 
###### 

################################################################## 


rep geos_das@dao : /product ion/geos_das/geosl . 1/gcm/ vc5 . 2b20 .macro vc5 . 2b20 

rep geos_das<9dao : /production/ geos_das/ geosl . 1/ applications/ gemprodn . f . 
rep geos_dasQdao:/production/geos_das/geosl . l/applications/progprs .f . 
rep geos_dasQdao : /product ion/geos_das/geosl . 1/applications/progsig. f . 
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rep geos_das<Odao:/production/geos_das/geosl . 1/applicat ions/diagprs .1 . 

rep geos_das©dao: /product ion/geos„das/geosl. i/applications/vrslvb20. com vrslv.com 

rep geos_das©dao: /product ion/geos_das/geosl . 1/applications/vstate.com . 

rep geos_das©dao:/production/geos_das/geosl. 1/applications/vreal . com . 

rep geos_das©dao: /production/geos_das/geosl . 1/applications/vphyscom. com . 

rep geos_das©dao: /production/geos_das/geosl . 1/applications/vmetric. com . 

rep geos_das©dao: /production/geos_das/geosl . 1/applications/metricc. com . 

rep geos_das©dao : /production/geos_das/geosl . 1/applicat ions/vbudget . com . 

rep geos_das©dao: /product ion/geos_das/geosl. 1/applicat ions/pblparm. com . 

rep geos_das©dao:/production/geos_das/geosl. 1/applicat ions/iau. com . 

rep geos_das©dao : /production/geos_das/geosl . 2/data/bcl44091 . yl985 gembes 
rep geos_das©dao : /product ion/geos_das/geosl . 2/data/gcmo3 . data . 
rep geos_das®dao : /product ion/geos_das/geosl .2/data/zonalT. data . 


12.4 Compile User Application and Output Routines 

Since the GEOS GCM library only contains subroutines from GCMDRV and below, the User 
must provide compiled Main Programs, or Applications, to perform specific experiments. 
For complete consistancy with the GEOS library, the Application and Output routines 
should be compiled in the same manner as the background model. An inline compile routine 
to accomplish this is shown below: 


################################################################## 

###### 

###### Create COMPILE Routine and Compile Userappl & Output 
###### 

################################################################## 

/bin/rm compile 
cat > compile CC’EOF’ 

/lib/epp -P -I . $1.1 $l.cpp 

fpp -dec $i.cpp > $l.fpp 

clt77 -i 64 -a static -V -exm $l.lpp 

rav $l.fpp.o $l.o 

'EOF* 

chraod +x compile 

./compile gemprodn 
./compile progprs 
./compile progs ig 
./compile diagprs 
./compile diagsize 


93 



Here we have compiled those routines which were remote copied in Section 12.3. 


12.5 Load and Run GEOS Simulation 


Following successful compilations, the User must load the Application and Output routines 
together with the GEOS GCM library, creating an executable a.out. Then, assuming that 
an appropriate GEOS GCM Restart has been brought to the current local directory as 
gcmrst , the GEOS simulation may be performed: 

################################################################## 

###### 

###### Load GCM, Applications and Output Routines, and Run 
###### 

################################################################## 


segldr -M,s -1 . /vc5.2b20 — gcmprodn.o progs ig.o progprs.o diagprs.o diagsize.o 


assign -a gcmrun. namelist 
assign -a gcmrst 


fort . 11 
fort . 12 


assign -a 
assign -a 
assign -a 
assign -a 
assign -a 

assign -a 
assign -a 
assign -a 


prog. sig 
prog.prs 
diag.prs.a 
diag.prs .b 
diag.prs.c 

gcmbcs 
gcmo3.dat a 
zonalT.data 


-N 

ieee 

-F 

f 77 

-N 

ieee 

-F 

f 77 

-N 

ieee 

-F 

f 77 

-N 

ieee 

-F 

f 77 

-N 

ieee 

-F 

f 77 

-N 

ibm 

-F 

COS 


fort . 34 
fort . 35 
fort .41 
fort .42 
fort .43 

fort . 16 
fort . 20 
fort . 52 


setenv NCPUS 4 
a.out 


Note that NCPUS has been set as an environmental variable equal to 4 in order to take 
advantage of the macro-tasking properties of the GEOS GCM. In this example, file units 34, 
35, 41, 42, 43 are associated with Output data streams, and file units 16, 20, 52 are input 
datasets and boundary conditions for the GEOS GCM. Once the GEOS GCM simulation 
has been completed, the resulting output data streams should be moved to permanent 
storage. 
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